Project overview

The goal of the Pyromania project is to test how terrestrial subsides (plant biomass loading) and burning influence DOM composition and degradation. We used a manipulative experiment to assess a range of plant material quantities (0-400g per tank) and fire treatment (burned vs unburned material) and the non-linearity of these effects on aquatic systems through 4 time-point samplings. We used 400L aquatic mesocosms and ran the experiment for ~90d in 2021-2022.

DATA SETS
This data set is among 3 to be generated for the project and focuses on:

  1. dissolved organic carbon (DOC) and
  2. Fluorescence and UV-vis spectroscopic absorbance indices
  3. DOC Degradation from either photo- or microbial decomposition
  4. plant biomass decomposition as % mass loss

TIME POINTS

  • Time-0 (T0) = before the addition of plant materials, plankton were stocked at this stage
  • Time-1 (T1) = 10 days after the addition of plant materials
  • Time-2 (T2) = 31 days after …
  • Time-3 (T3) = 59 days after …
  • Time-4 (T4) = 89 days after …

###DOC
import and for loop to clean up all files and make stacked data in single df

## import treatment IDs
IDs<-read.csv("data/treatment.IDs.csv")

##### grab files in a list
Total.DOC.files <- list.files(path="data/DOC.TN", pattern = "csv$", full.names = T)
Total.DOC.files
## [1] "data/DOC.TN/DOC_T0.csv" "data/DOC.TN/DOC_T1.csv" "data/DOC.TN/DOC_T2.csv"
## [4] "data/DOC.TN/DOC_T3.csv" "data/DOC.TN/DOC_T4.csv"
##### what are the file names, sans extensions using package 'tools'
file.names<-file_path_sans_ext(list.files(path="data/DOC.TN", pattern = "csv$", full.names = F))
file.names
## [1] "DOC_T0" "DOC_T1" "DOC_T2" "DOC_T3" "DOC_T4"
############ formatting all data in for loop
  for(i in 1:length(Total.DOC.files))
    {
  data<-read.csv(Total.DOC.files[i], sep=",")
  data<-data[,c(1,3,4)] # removed columns we don't need
  data$File<-Total.DOC.files[i]
  colnames(data)<-c("Tank", "DOC..mg.L", "TN..mg.L", "File")
  data$Tank<- IDs$Tank
  data$Tank<-as.numeric(as.character(data$Tank)) # make the column of samples all numeric
  data <- data[!is.na(as.numeric(as.character(data$Tank))),] # remove all rows that aren't numeric/tanks
  data$Treatment<-IDs$Treatment
  data$plant.mass..g<-IDs$plant.mass..g
  make.names(assign(paste(file.names[i], sep=""), data)) # make the file name the name of new df for loop df
  }
########## this is the end of the loop


ls() #see all dfs you've made, the above will be df matching their file names
##  [1] "data"            "DOC_T0"          "DOC_T1"          "DOC_T2"         
##  [5] "DOC_T3"          "DOC_T4"          "Fig.formatting"  "file.names"     
##  [9] "i"               "IDs"             "Total.DOC.files"
#Combine files from loop to single df
DOC.df<-rbind(DOC_T0, DOC_T1, DOC_T2, DOC_T3, DOC_T4)

DOC.df$File <- sapply(strsplit(DOC.df$File, "/"), `[`, 3) # extract sample names

# alternative way to code the above
#give the 10th-24th character of the file name, removing the rest
#DOC.df$File<-substr(DOC.df$File, 13, 27) 

#alternatively
# remove the 9 letters ('^.) at start 
# remove the 4 letters (.$') at end
#DOC.df$File<-gsub('^.........|....$', '', DOC.df$File) 

# if equals DOC_T0_11052021 then, T0, if not then T1
DOC.df$Time.point<- as.factor(ifelse(DOC.df$File=="DOC_T0.csv", "T0",
         ifelse (DOC.df$File=="DOC_T1.csv", "T1",
           ifelse (DOC.df$File=="DOC_T2.csv", "T2", 
                   ifelse(DOC.df$File=="DOC_T3.csv", "T3", "T4")))))

#rearrange
DOC.df<- DOC.df %>% 
  select(File, Time.point, Treatment, Tank, plant.mass..g, DOC..mg.L, TN..mg.L) 

DOC.df$Treatment<-as.factor(DOC.df$Treatment)

DOC analysis:

  • notes on GAMS…. Note that fitting via restricted maximum likelihood (REML) is advised to give stable results. Fitting the smoothing parameter (sp) to user defined levels not advisable (but certainly is done). Number of basis functions important as well, small # of basis functions (k) will limit wiggliness (and make the model more linear than should be). This wiggliness is set by using the k values. Don’t want too high, or too low…

When using the non-linear smoothing, this is the s(xxxx). Can add a new predictor by adding another smooth function with s(x1) + s(x2). When the variable is inside the smooth function, this is to account for the nonlinear shape. Two smoothing together and you have ADDITIVE Nonlinear modeling. Not all variables need to be non-linear, but adding a predictor outside of the s() then this is allowing for linearity.

In practice, continuous variables are rarely linear in GAMs, automatic smoothing will make a linear shape if setting high smoothing function (ie.. s =1000), then a non-linear variable appears linear. However, setting categorical variables as predicots linear is more common.

Factor-smooth interaction also an option, this is when s(x1 by = fac), this sets different smoothing for different variables of “fac”. Usually, the different smoothers are combined with a varying intercept incase the different categories are different in means and slopes, this would be by adding the s(x1 by = fac) + fac, wheere the +fac allows for a different slope

EDF - effective degrees of freedom, equate with wiggliness, with edf =1 as a straight line, and higher edfs as more wiggly. GAM smoother significance is described as not being able to draw a horizontal line through the data. Always check concurvity, which is the collinearity with models from 0-1.

######## T0 model
m1.DOC.T0.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")

m2.DOC.T0.gam=gam(DOC..mg.L ~  Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")

m3.DOC.T0.gam=gam(DOC..mg.L ~  s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")

T0.DOC.AIC<-AIC(m1.DOC.T0.gam, m2.DOC.T0.gam, m3.DOC.T0.gam)

DOC.sum_T0=summary(m3.DOC.T0.gam)
DOC.sum_T0=anova.gam(m3.DOC.T0.gam)
gam.check(m3.DOC.T0.gam, rep=1000)
draw(m3.DOC.T0.gam)
concrvity(m3.DOC.T0.gam)
par(mfrow = c(2, 2))
plot(m3.DOC.T0.gam, all.terms = TRUE, page=1)


DOC.diff.T0<-plot_difference(
  m1.DOC.T0.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
) +
  coord_cartesian(ylim = c(-10,13))+
  Fig.formatting

###########  
#plot for the model output on rawdata
DOC.T0.mod.plot<-
  plot_smooths(
  model = m3.DOC.T0.gam,
  series = plant.mass..g
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),], 
             aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  coord_cartesian(ylim=c(0, 60)) +
  ggtitle("Day-0") +
  ylab("DOC (mg/L)") +
  xlab("plant material (g)") +
  Fig.formatting

######## T1 model
m1.DOC.T1.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")

m2.DOC.T1.gam=gam(DOC..mg.L ~  Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")

m3.DOC.T1.gam=gam(DOC..mg.L ~  s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")

T1.DOC.AIC<-AIC(m1.DOC.T1.gam, m2.DOC.T1.gam, m3.DOC.T1.gam)

summary(m1.DOC.T1.gam)
T1.DOC.anova=anova(m1.DOC.T1.gam)
gam.check(m1.DOC.T1.gam, rep=1000)
draw(m1.DOC.T1.gam)
concrvity(m1.DOC.T1.gam)
par(mfrow = c(2, 2))
plot(m1.DOC.T1.gam, all.terms = TRUE, page=1)


# model predictions
DOC.diff.T1<-plot_difference(
  m1.DOC.T1.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-10,13))+
  Fig.formatting

###########  
#plot for the model output on rawdata
DOC.T1.mod.plot<-
  plot_smooths(
  model = m1.DOC.T1.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  coord_cartesian(ylim=c(0, 60)) +
  ggtitle("Day-10") +
  ylab("DOC (mg/L)") +
  xlab("plant material (g)") +
  Fig.formatting

########## T2
m1.DOC.T2.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")

m2.DOC.T2.gam=gam(DOC..mg.L ~  Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")

m3.DOC.T2.gam=gam(DOC..mg.L ~  s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")

T2.DOC.AIC<-AIC(m1.DOC.T2.gam, m2.DOC.T2.gam, m3.DOC.T2.gam)

summary(m1.DOC.T2.gam)
T2.DOC.anova=anova.gam(m1.DOC.T2.gam)
gam.check(m1.DOC.T2.gam, rep=1000)
draw(m1.DOC.T2.gam)
concrvity(m1.DOC.T2.gam)
par(mfrow = c(2, 2))
plot(m1.DOC.T2.gam, all.terms = TRUE, page=1)

# model predictions
DOC.diff.T2<-plot_difference(
  m1.DOC.T2.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-10,13))+
  Fig.formatting

###########  
#plot for the model output on rawdata
DOC.T2.mod.plot<-
  plot_smooths(
  model = m1.DOC.T2.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  coord_cartesian(ylim=c(0, 60)) +
  ggtitle("Day-31") +
  ylab("DOC (mg/L)") +
  xlab("plant material (g)") +
  Fig.formatting

########## T3
m1.DOC.T3.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")

m2.DOC.T3.gam=gam(DOC..mg.L ~  Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")

m3.DOC.T3.gam=gam(DOC..mg.L ~  s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")

T3.DOC.AIC<-AIC(m1.DOC.T3.gam, m2.DOC.T3.gam, m3.DOC.T3.gam)

summary(m1.DOC.T3.gam)
T3.DOC.anova=anova.gam(m1.DOC.T3.gam)
gam.check(m1.DOC.T3.gam, rep=1000)
draw(m1.DOC.T3.gam)
concrvity(m1.DOC.T3.gam)
par(mfrow = c(2, 2))
plot(m1.DOC.T3.gam, all.terms = TRUE, page=1)

# model predictions
DOC.diff.T3<-plot_difference(
  m1.DOC.T3.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-10,13))+
  Fig.formatting

###########  
#plot for the model output on rawdata
DOC.T3.mod.plot<-
  plot_smooths(
  model = m1.DOC.T3.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),], 
             aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  coord_cartesian(ylim=c(0, 60)) +
  ggtitle("Day-59") +
  ylab("DOC (mg/L)") +
  xlab("plant material (g)") +
  Fig.formatting

########## T4
m1.DOC.T4.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")

m2.DOC.T4.gam=gam(DOC..mg.L ~  Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")

m3.DOC.T4.gam=gam(DOC..mg.L ~  s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")

T4.DOC.AIC<-AIC(m1.DOC.T4.gam, m2.DOC.T4.gam, m3.DOC.T4.gam)

summary(m3.DOC.T4.gam)
T4.DOC.anova=anova(m1.DOC.T4.gam)
gam.check(m3.DOC.T4.gam, rep=1000)
draw(m3.DOC.T4.gam)
concrvity(m3.DOC.T4.gam)
par(mfrow = c(2, 2))
plot(m3.DOC.T4.gam, all.terms = TRUE, page=1)

###########  
#plot for the model output on rawdata
DOC.T4.mod.plot<-
  plot_smooths(
  model = m3.DOC.T4.gam,
  series = plant.mass..g
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),], 
             aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  coord_cartesian(ylim=c(0, 60)) +
  ggtitle("Day-89") +
  ylab("DOC (mg/L)") +
  xlab("plant material (g)") +
  Fig.formatting

# model predictions
DOC.diff.T4<-plot_difference(
  m1.DOC.T4.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-10,13))+
  Fig.formatting

#labels
mod.DOC.rep<-rep(c(
              "Day 10: Treatment with by-factor smooth", 
              "Day 10: Treatment with global smooth", 
              "Day 10: global smooth",
              "Day 31: Treatment with by-factor smooth", 
              "Day 31: Treatment with global smooth", 
              "Day 31: global smooth",
              "Day 59: Treatment with by-factor smooth", 
              "Day 59: Treatment with global smooth", 
              "Day 59: global smooth",
              "Day 89: Treatment with by-factor smooth", 
              "Day 89: Treatment with global smooth", 
              "Day 89: global smooth"))

#adding columns
mod.DOC.column<- data.frame(mod.DOC.rep)
AIC.DOC<-bind_rows(T1.DOC.AIC, T2.DOC.AIC, T3.DOC.AIC, T4.DOC.AIC)
AIC.DOC.mod<-as.data.frame(cbind(mod.DOC.column, AIC.DOC))

names(AIC.DOC.mod)[1]="Model"

AIC.DOC.mod <- AIC.DOC.mod %>%
  # Creating an empty column:
  add_column(Variable= c("DOC concentration"), .before="Model")


write.csv(AIC.DOC.mod, "output/AIC models/AIC.DOC.csv")

#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/DOC.mod.AIC.pdf", width = 6.8, height = 3.8)       # Export PDF
grid.table(AIC.DOC.mod, rows = NULL)
dev.off()

Generate dataframe for model outputs

####Parametric output
DOC.para.sums=as.data.frame(rbind(T1.DOC.anova$pTerms.table, T2.DOC.anova$pTerms.table, T3.DOC.anova$pTerms.table, T4.DOC.anova$pTerms.table))
names(DOC.para.sums)[1]="df/edf"

#create empty column for ref.df
DOC.para.sums <- DOC.para.sums %>%
  # Creating an empty column:
  add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
  add_column(Ref.df = "-", .after="df/edf")
  
#####Smooth output
DOC.smooth.sums=as.data.frame(rbind(T1.DOC.anova$s.table, T2.DOC.anova$s.table, T3.DOC.anova$s.table, T4.DOC.anova$s.table))
names(DOC.smooth.sums)[1]="df/edf"
#create empty column for ref.df
DOC.smooth.sums <- DOC.smooth.sums %>%
  # Creating an empty column:
  add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")

#merge the dataframes
DOC.mod.allsums=rbind(DOC.para.sums,DOC.smooth.sums)
rownames(DOC.mod.allsums)<-c(1:12)

#create new ID column
DOC.mod.allsums <- DOC.mod.allsums %>%
  # Creating an empty column:
  add_column(Variable= c("DOC concentration"), .before="Model")

#convert dataframe into a table for export
pdf("output/Codys.tables/DOC.mod.allsums.pdf", width = 9.3, height = 3.8)       # Export PDF
grid.table(DOC.mod.allsums, rows= NULL)
dev.off()
## quartz_off_screen 
##                 2

DOC: compile raw plots and model-diff plots for final figures

###### compile the plots with effect plots}

extract.legend <- get_legend(
  # create some space to the left of the legend
  DOC.T1.mod.plot + theme(legend.box.margin = margin(0, 0, 0, 10)))


DOC.alltimes<-plot_grid(
  DOC.T0.mod.plot+ theme(legend.position = "none"),
  DOC.T1.mod.plot+ theme(legend.position = "none"),
  DOC.T2.mod.plot+ theme(legend.position = "none"),
  DOC.T3.mod.plot+ theme(legend.position = "none"),
  DOC.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
  rel_widths = c(8,8,8,8,8,3), ncol=6)

ggsave("figures/DOC.alltimes.mod.pdf", height=4, width=18)
DOC.alltimes
**Figure**. Dissolved organic carbon (DOC) concentration across time in treatments receiving burned and unburned plant material at the start of the experiment and four sampling periods after plant material added. Black lines with gray confidence intervals indicate global smoothers across all data points; solid (burned) and dotted (unburned) black lines together represent treatment-level intercepts with global smoothers; colored lines indicate factor-smooths that vary between treatments.

Figure. Dissolved organic carbon (DOC) concentration across time in treatments receiving burned and unburned plant material at the start of the experiment and four sampling periods after plant material added. Black lines with gray confidence intervals indicate global smoothers across all data points; solid (burned) and dotted (unburned) black lines together represent treatment-level intercepts with global smoothers; colored lines indicate factor-smooths that vary between treatments.

DOC.mod.diffs<-plot_grid(
  DOC.diff.T0+ theme(legend.position = "none")+ ggtitle("Day-0"),
  DOC.diff.T1+ theme(legend.position = "none")+ ggtitle("Day-10"),
  DOC.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
  DOC.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
  DOC.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
  rel_widths = c(8,8,8,8,8), ncol=5)
DOC.mod.diffs
**Figure**. Model effects from GAMs with differences between smoothers for DOC concentration across time in treatments receiving burned and unburned plant material. The areas in pink show where there are significant differences between the two smoothers, indicating treatment effects.

Figure. Model effects from GAMs with differences between smoothers for DOC concentration across time in treatments receiving burned and unburned plant material. The areas in pink show where there are significant differences between the two smoothers, indicating treatment effects.

ggsave("figures/DOC.mod.diffs.pdf", height=4, width=16)

UVBIO files: DOC degradation

##### grab files in a list
UVBio.files <- list.files(path="/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO", pattern = "csv$", full.names = T)
UVBio.files
## [1] "/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO/DOC_NoUV_T1.csv"
## [2] "/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO/DOC_NoUV_T2.csv"
## [3] "/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO/DOC_NoUV_T3.csv"
## [4] "/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO/DOC_UV_T1.csv"  
## [5] "/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO/DOC_UV_T2.csv"  
## [6] "/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO/DOC_UV_T3.csv"
##### what are the file names, sans extensions using package 'tools'
file.names<-file_path_sans_ext(list.files(path="/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO", pattern = "csv$", full.names = F))
file.names
## [1] "DOC_NoUV_T1" "DOC_NoUV_T2" "DOC_NoUV_T3" "DOC_UV_T1"   "DOC_UV_T2"  
## [6] "DOC_UV_T3"
############ formatting all data in for loop
  for(i in 1:length(UVBio.files))
    {
  data<-read.csv(UVBio.files[i], sep=",") #skip 13 rows where standards are
  data<-data[,c(1,3,4)] # removed columns we don't need
  data$File<-UVBio.files[i]
  colnames(data)<-c("Tank", "DOC..mg.L", "TN..mg.L", "File") 
  data$Tank<-as.character(data$Tank)
  data$Treatment<-IDs$Treatment # add treatment info
  data$plant.mass..g<-IDs$plant.mass..g # add plant biomass info
  make.names(assign(paste(file.names[i], sep=""), data)) # make the file name the name of new df for loop df
  }
########## this is the end of the loop

ls() #see all dfs you've made, the above will be df matching their file names
##  [1] "AIC.DOC"         "AIC.DOC.mod"     "data"            "DOC_NoUV_T1"    
##  [5] "DOC_NoUV_T2"     "DOC_NoUV_T3"     "DOC_T0"          "DOC_T1"         
##  [9] "DOC_T2"          "DOC_T3"          "DOC_T4"          "DOC_UV_T1"      
## [13] "DOC_UV_T2"       "DOC_UV_T3"       "DOC.alltimes"    "DOC.df"         
## [17] "DOC.diff.T0"     "DOC.diff.T1"     "DOC.diff.T2"     "DOC.diff.T3"    
## [21] "DOC.diff.T4"     "DOC.mod.allsums" "DOC.mod.diffs"   "DOC.para.sums"  
## [25] "DOC.smooth.sums" "DOC.sum_T0"      "DOC.T0.mod.plot" "DOC.T1.mod.plot"
## [29] "DOC.T2.mod.plot" "DOC.T3.mod.plot" "DOC.T4.mod.plot" "extract.legend" 
## [33] "Fig.formatting"  "file.names"      "i"               "IDs"            
## [37] "m1.DOC.T0.gam"   "m1.DOC.T1.gam"   "m1.DOC.T2.gam"   "m1.DOC.T3.gam"  
## [41] "m1.DOC.T4.gam"   "m2.DOC.T0.gam"   "m2.DOC.T1.gam"   "m2.DOC.T2.gam"  
## [45] "m2.DOC.T3.gam"   "m2.DOC.T4.gam"   "m3.DOC.T0.gam"   "m3.DOC.T1.gam"  
## [49] "m3.DOC.T2.gam"   "m3.DOC.T3.gam"   "m3.DOC.T4.gam"   "mod.DOC.column" 
## [53] "mod.DOC.rep"     "T0.DOC.AIC"      "T1.DOC.AIC"      "T1.DOC.anova"   
## [57] "T2.DOC.AIC"      "T2.DOC.anova"    "T3.DOC.AIC"      "T3.DOC.anova"   
## [61] "T4.DOC.AIC"      "T4.DOC.anova"    "Total.DOC.files" "UVBio.files"
#Combine files from loop to single df
UVBIO.df<-rbind(DOC_NoUV_T1, DOC_UV_T1, DOC_NoUV_T2, DOC_UV_T2, DOC_NoUV_T3, DOC_UV_T3)

# make new factors for filtered (F/NF) and UV (+/-)
# ** caution! Make sure this matches sheet layout after the loop!

UVBIO.df$Filter.Trt<-rep(c("Filtered", "Unfiltered"), each=30) # repeat 30 times
UVBIO.df$UV.Trt<-rep(c("NoUV", "UV"), each=60) # repeat 60 times

# check tank names in DF and make sure they are what you think they are (i.e 1-30 but w/ weird names)
Tank.fix<-rep(c(1:30), times=12) # repeat 12 times to fix all the strange tank names 

#compare tank IDS
df.test<-cbind(Tank.fix, UVBIO.df$Tank) # yep!

UVBIO.df$Tank<-as.factor(Tank.fix) #replace old tank names

# extract out the file names
#Extract file
UVBIO.df$File=sapply(strsplit(UVBIO.df$File, "/"), `[`, 9)
#take out redundancies
UVBIO.df$DOC..mg.L<-gsub("NPOC:","", UVBIO.df$DOC..mg.L)
UVBIO.df$DOC..mg.L<- gsub("mg/L", "", UVBIO.df$DOC..mg.L)
UVBIO.df$TN..mg.L<-gsub("TN:","", UVBIO.df$TN..mg.L)
UVBIO.df$TN..mg.L<- gsub("mg/L", "", UVBIO.df$TN..mg.L)

#convert DOC and TN to numeric
UVBIO.df$DOC..mg.L<-as.numeric(as.character(UVBIO.df$DOC..mg.L))
UVBIO.df$TN..mg.L<- as.numeric(as.character(UVBIO.df$TN..mg.L))

# add a time point
# if equals DOC_UV_T1_111721 then, T1, if not then "nothing"
UVBIO.df$Time.point<- as.factor(ifelse(UVBIO.df$File=="DOC_UV_T1.csv", "T1",
  ifelse(UVBIO.df$File=="DOC_NoUV_T1.csv", "T1",
          ifelse(UVBIO.df$File=="DOC_UV_T2.csv", "T2",
                ifelse(UVBIO.df$File=="DOC_NoUV_T2.csv", "T2", "T3")))))

#rearrange
UVBIO.df<- UVBIO.df %>% 
  dplyr::select(File, Time.point, Treatment, Filter.Trt, UV.Trt, Tank, plant.mass..g, DOC..mg.L, TN..mg.L)


# make levels for plotting and stats
# add in starting DOC and TN
# original data is... DOC_T0, DOC_T1, DOC_T2, DOC_T3, DOC_T4 -- but no need for DOC T0 and T4
# repeat 4x to account for the 4 treatments in each tank
DOC_T1.start<-rep(DOC_T1$DOC..mg.L, times=4)
DOC_T2.start<-rep(DOC_T2$DOC..mg.L, times=4)
DOC_T3.start<-rep(DOC_T3$DOC..mg.L, times=4)

TN_T1.start<-rep(DOC_T1$TN..mg.L, times=4)
TN_T2.start<-rep(DOC_T2$TN..mg.L, times=4)
TN_T3.start<-rep(DOC_T3$TN..mg.L, times=4)

# compile and add in as a new column
UVBIO.df$Start.DOC<-c(DOC_T1.start, DOC_T2.start, DOC_T3.start)
UVBIO.df$Start.TN<-c(TN_T1.start, TN_T2.start, TN_T3.start)

UVBIO.df$Burn.Filt.UV<-interaction(UVBIO.df$Treatment, UVBIO.df$Filter.Trt, UVBIO.df$UV.Trt)
UVBIO.df$Filt.UV<-interaction(UVBIO.df$Filter.Trt, UVBIO.df$UV.Trt)

###SUVA SUVA calculations Time point 0

#load file and rearrange

suva.t0.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/Results_T0.csv")
suva.t0.df$Sample.Name=gsub("\\..*", "", suva.t0.df$Sample.Name) 
names(suva.t0.df)[1]="Tank"
suva.t0.df=suva.t0.df[-1,]
suva.t0.df$Tank=gsub("T", "", suva.t0.df$Tank) 
suva.t0.df$Tank=as.numeric(suva.t0.df$Tank)
suva.t0.df=suva.t0.df[order(suva.t0.df$Tank),]

#add important info
suva.t0.df$DOC..mg.L=DOC_T0$DOC..mg.L
suva.t0.df$Time.point="T0"
suva.t0.df$Treatment=IDs$Treatment
suva.t0.df$plant.mass..g=IDs$plant.mass..g

#calculate SUVA
suva.t0.df$SUVA=(suva.t0.df[suva.t0.df$Tank,]$abs254/suva.t0.df[suva.t0.df$Tank,]$DOC..mg.L)*100

Time point 1

#load file and rearrange
library(dplyr)

suva.t1.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/Results_T1.csv")
suva.t1.df=suva.t1.df[-1,]
suva.t1.df=suva.t1.df[-1,]
names(suva.t1.df)[1]="Tank"
suva.t1.df$Tank=gsub("TA", "", suva.t1.df$Tank) 
suva.t1.df$Tank=sub("(.{2})(.*)", "\\1.\\2", suva.t1.df$Tank)
suva.t1.df$Tank=sub("\\..*", "", suva.t1.df$Tank)
suva.t1.df$Tank=as.numeric(suva.t1.df$Tank)
suva.t1.df=suva.t1.df[order(suva.t1.df$Tank),]

#add important info
suva.t1.df$DOC..mg.L=DOC_T1$DOC..mg.L
suva.t1.df$Time.point="T1"
suva.t1.df$Treatment=IDs$Treatment
suva.t1.df$plant.mass..g=IDs$plant.mass..g

#calculate SUVA
suva.t1.df$SUVA=(suva.t1.df[suva.t1.df$Tank,]$abs254/suva.t1.df[suva.t1.df$Tank,]$DOC..mg.L)*100

Time point 2

#load file and rearrange
suva.t2.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/Results_T2_d1A.edit.csv")
names(suva.t2.df)[1]="Tank"
suva.t2.df$Tank=gsub("TA", "", suva.t2.df$Tank) 
suva.t2.df$Tank=sub("(.{2})(.*)", "\\1.\\2", suva.t2.df$Tank)
suva.t2.df$Tank=sub("\\..*", "", suva.t2.df$Tank)
suva.t2.df$Tank=as.numeric(suva.t2.df$Tank)
suva.t2.df=suva.t2.df[order(suva.t2.df$Tank),]

#add important info
suva.t2.df$DOC..mg.L=DOC_T2$DOC..mg.L
suva.t2.df$Time.point="T2"
suva.t2.df$Treatment=IDs$Treatment
suva.t2.df$plant.mass..g=IDs$plant.mass..g

#calculate SUVA
suva.t2.df$SUVA=(suva.t2.df[suva.t2.df$Tank,]$abs254/suva.t2.df[suva.t2.df$Tank,]$DOC..mg.L)*100

Time point 3

suva.t3.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/Results_T3.csv")
names(suva.t3.df)[1]="Tank"
suva.t3.df$Tank=gsub("TA", "", suva.t3.df$Tank) 
suva.t3.df$Tank=sub("(.{2})(.*)", "\\1.\\2", suva.t3.df$Tank)
suva.t3.df$Tank=sub("\\..*", "", suva.t3.df$Tank)
suva.t3.df$Tank=as.numeric(suva.t3.df$Tank)
suva.t3.df=suva.t3.df[order(suva.t3.df$Tank),]


#add important info
suva.t3.df$DOC..mg.L=DOC_T3$DOC..mg.L
suva.t3.df$Time.point="T3"
suva.t3.df$Treatment=IDs$Treatment
suva.t3.df$plant.mass..g=IDs$plant.mass..g
class(suva.t3.df$DOC..mg.L)
## [1] "numeric"
#calculate SUVA
suva.t3.df$SUVA=(suva.t3.df[suva.t3.df$Tank,]$abs254/suva.t3.df[suva.t3.df$Tank,]$DOC..mg.L)*100

Time point 4

suva.t4.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/Results_T4.csv")
names(suva.t4.df)[1]="Tank"
suva.t4.df$Tank=gsub("TA", "", suva.t4.df$Tank) 
suva.t4.df$Tank=sub("(.{2})(.*)", "\\1.\\2", suva.t4.df$Tank)
suva.t4.df$Tank=sub("\\..*", "", suva.t4.df$Tank)
suva.t4.df$Tank=as.numeric(suva.t4.df$Tank)
suva.t4.df=suva.t4.df[order(suva.t4.df$Tank),]

#add important info
suva.t4.df$DOC..mg.L=DOC_T4$DOC..mg.L
suva.t4.df$Time.point="T4"
suva.t4.df$Treatment=IDs$Treatment
suva.t4.df$plant.mass..g=IDs$plant.mass..g

#calculate SUVA
suva.t4.df$SUVA=(suva.t4.df[suva.t4.df$Tank,]$abs254/suva.t4.df[suva.t4.df$Tank,]$DOC..mg.L)*100

Combine all SUVA time points

suva.df=rbind(suva.t0.df, suva.t1.df, suva.t2.df, suva.t3.df, suva.t4.df)

Other indices

# Time 0
index.T0.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/T0_indices.csv")
index.T0.df$Time.point="T0"
names(index.T0.df)[6]="Humification.Index"
index.T0.df=index.T0.df[-c(7:11,13,14)]

#Time 1
index.T1.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/T1_indices.csv")
index.T1.df$Time.point="T1"
index.T1.df=index.T1.df[-c(7:11,13,14)]

#Time 2
index.T2.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/T2_indices.csv")
index.T2.df$Time.point="T2"
index.T2.df=index.T2.df[-c(7:11,13,14)]

#Time 3
index.T3.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/T3_indices.csv")
index.T3.df$Time.point="T3"
index.T3.df=index.T3.df[-c(7:11,13,14)]


#Time 4
index.T4.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/T4_indices.csv")
index.T4.df$Time.point="T4"
index.T4.df=index.T4.df[-c(7:11,13,14)]

#combine all index time points
index.df=rbind(index.T0.df, index.T1.df, index.T2.df, index.T3.df, index.T4.df)

Add these new metrics to DOC df import and do a loop to clean up all files and make stacked data in single df

#Combine files from loop to single df
DOC_T0$SUVA=suva.df[suva.df$Time.point=="T0",]$SUVA
DOC_T1$SUVA=suva.df[suva.df$Time.point=="T1",]$SUVA
DOC_T2$SUVA=suva.df[suva.df$Time.point=="T2",]$SUVA
DOC_T3$SUVA=suva.df[suva.df$Time.point=="T3",]$SUVA
DOC_T4$SUVA=suva.df[suva.df$Time.point=="T4",]$SUVA

#HIX
DOC_T0$hix=index.df[index.df$Time.point=="T0",]$Humification.Index
DOC_T1$hix=index.df[index.df$Time.point=="T1",]$Humification.Index
DOC_T2$hix=index.df[index.df$Time.point=="T2",]$Humification.Index
DOC_T3$hix=index.df[index.df$Time.point=="T3",]$Humification.Index
DOC_T4$hix=index.df[index.df$Time.point=="T4",]$Humification.Index

#Fluorescence
DOC_T0$fl=index.df[index.df$Time.point=="T0",]$Fluroscence.Index
DOC_T1$fl=index.df[index.df$Time.point=="T1",]$Fluroscence.Index
DOC_T2$fl=index.df[index.df$Time.point=="T2",]$Fluroscence.Index
DOC_T3$fl=index.df[index.df$Time.point=="T3",]$Fluroscence.Index
DOC_T4$fl=index.df[index.df$Time.point=="T4",]$Fluroscence.Index

#slope ratio
DOC_T0$sr=suva.df[suva.df$Time.point=="T0",]$Spectral.Slope.Ratio
DOC_T1$sr=suva.df[suva.df$Time.point=="T1",]$Spectral.Slope.Ratio
DOC_T2$sr=suva.df[suva.df$Time.point=="T2",]$Spectral.Slope.Ratio
DOC_T3$sr=suva.df[suva.df$Time.point=="T3",]$Spectral.Slope.Ratio
DOC_T4$sr=suva.df[suva.df$Time.point=="T4",]$Spectral.Slope.Ratio

#freshness (BIX)
DOC_T0$fr=index.df[index.df$Time.point=="T0",]$Freshness.Index
DOC_T1$fr=index.df[index.df$Time.point=="T1",]$Freshness.Index
DOC_T2$fr=index.df[index.df$Time.point=="T2",]$Freshness.Index
DOC_T3$fr=index.df[index.df$Time.point=="T3",]$Freshness.Index
DOC_T4$fr=index.df[index.df$Time.point=="T4",]$Freshness.Index

#combine
DOC.df<-rbind(DOC_T0, DOC_T1, DOC_T2, DOC_T3, DOC_T4)

#Extract file
DOC.df$File <- sapply(strsplit(DOC.df$File, "/"), `[`, 3) # extract sample names

# alternative way to code the above
#give the 10th-24th character of the file name, removing the rest
#DOC.df$File<-substr(DOC.df$File, 13, 27) 

#alternatively
# remove the 9 letters ('^.) at start 
# remove the 4 letters (.$') at end
#DOC.df$File<-gsub('^.........|....$', '', DOC.df$File) 

# if equals DOC_T0 then, T0, if not then T1..etc
DOC.df$Time.point<- as.factor(ifelse(DOC.df$File=="DOC_T0.csv", "T0",
         ifelse (DOC.df$File=="DOC_T1.csv", "T1",
           ifelse (DOC.df$File=="DOC_T2.csv", "T2", 
                   ifelse(DOC.df$File=="DOC_T3.csv", "T3", "T4")))))

#rearrange
library(dplyr)
DOC.df<- DOC.df %>% 
  dplyr::select(File, Time.point, Treatment, Tank, plant.mass..g, DOC..mg.L, TN..mg.L, SUVA, sr, hix, fl, fr) 

DOC.df=DOC.df%>%
    mutate_if(is.character,as.factor)

SUVA analysis and plots

######## T0 model
m1.SUVA.T0.gam=gam(SUVA ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")

m2.SUVA.T0.gam=gam(SUVA ~  Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")

m3.SUVA.T0.gam=gam(SUVA ~  s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")

T0.suva.AIC<-AIC(m1.SUVA.T0.gam, m2.SUVA.T0.gam, m3.SUVA.T0.gam)
T0.suva.AIC
##                      df      AIC
## m1.SUVA.T0.gam 8.544088 68.24603
## m2.SUVA.T0.gam 4.730649 70.94187
## m3.SUVA.T0.gam 3.000364 81.69088
summary(m1.SUVA.T0.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## SUVA ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         2.3007     0.1664  13.824 1.04e-12 ***
## Treatmentunburned   1.0159     0.2354   4.316 0.000251 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   1.000  1.000 0.995   0.329
## s(plant.mass..g):Treatmentunburned 3.734  4.544 1.598   0.176
## 
## R-sq.(adj) =   0.46   Deviance explained = 56.7%
## -REML = 34.112  Scale est. = 0.41545   n = 30
anova.gam(m1.SUVA.T0.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## SUVA ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F  p-value
## Treatment  1 18.63 0.000251
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   1.000  1.000 0.995   0.329
## s(plant.mass..g):Treatmentunburned 3.734  4.544 1.598   0.176
gam.check(m1.SUVA.T0.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-1.233511e-05,6.231324e-06]
## (score 34.11213 & scale 0.4154513).
## Hessian positive definite, eigenvalue range [1.233492e-05,13.1468].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 1.00    1.06    0.49
## s(plant.mass..g):Treatmentunburned 9.00 3.73    1.06    0.56
draw(m1.SUVA.T0.gam)

concrvity(m1.SUVA.T0.gam)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5   e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     4.60e-25
## 3 worst    s(plant.mass..g):Treatmentunburned   4.58e-25
## 4 observed para                                 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned     5.91e-30
## 6 observed s(plant.mass..g):Treatmentunburned   1.91e-29
## 7 estimate para                                 5   e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned   1.24e-27
par(mfrow = c(2, 2))
plot(m1.SUVA.T0.gam, all.terms = TRUE, page=1)

# model predictions
suva.diff.T0<-plot_difference(
  m1.SUVA.T0.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
suva.T0.mod.plot<-
  plot_smooths(
  model = m3.SUVA.T0.gam,
  series = plant.mass..g,
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),], 
             aes(x=plant.mass..g, y=SUVA, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  ggtitle("Day-0") +
  ylab(expression('SUVA'["254"]~'(mgC'~ L^-1~m^-1*')'))+
  coord_cartesian(ylim = c(0,12))+
  xlab("plant material (g)") +
  Fig.formatting

suva.T0.mod.plot

######## T1 model
m1.suva.T1.gam=gam(SUVA ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")

m2.suva.T1.gam=gam(SUVA ~  Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")

m3.suva.T1.gam=gam(SUVA ~  s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")

T1.suva.AIC<-AIC(m1.suva.T1.gam, m2.suva.T1.gam, m3.suva.T1.gam)
T1.suva.AIC
##                       df      AIC
## m1.suva.T1.gam 11.829790 16.68484
## m2.suva.T1.gam  7.758699 18.60925
## m3.suva.T1.gam  6.762722 17.68119
summary(m1.suva.T1.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## SUVA ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        2.28324    0.06804  33.559   <2e-16 ***
## Treatmentunburned  0.09731    0.09622   1.011    0.324    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F  p-value    
## s(plant.mass..g):Treatmentburned   4.597  5.541 20.70  < 2e-16 ***
## s(plant.mass..g):Treatmentunburned 3.354  4.095 14.88 6.98e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.857   Deviance explained = 90.1%
## -REML = 15.005  Scale est. = 0.069436  n = 30
T1.suva.anova=anova(m1.suva.T1.gam)
gam.check(m1.suva.T1.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [-6.613102e-06,3.156923e-06]
## (score 15.00512 & scale 0.06943626).
## Hessian positive definite, eigenvalue range [0.9040944,13.38198].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value  
## s(plant.mass..g):Treatmentburned   9.00 4.60    0.77   0.025 *
## s(plant.mass..g):Treatmentunburned 9.00 3.35    0.77   0.025 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
draw(m1.suva.T1.gam)

concrvity(m1.suva.T1.gam)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5   e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     4.60e-25
## 3 worst    s(plant.mass..g):Treatmentunburned   4.58e-25
## 4 observed para                                 5   e- 1
## 5 observed s(plant.mass..g):Treatmentburned     1.18e-29
## 6 observed s(plant.mass..g):Treatmentunburned   2.90e-30
## 7 estimate para                                 5   e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned   1.24e-27
par(mfrow = c(2, 2))
plot(m1.suva.T1.gam, all.terms = TRUE, page=1)

# model predictions
suva.diff.T1<-plot_difference(
  m1.suva.T1.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-5,5))+
  Fig.formatting

###########  
#plot for the model output on rawdata
suva.T1.mod.plot<-
  plot_smooths(
  model = m1.suva.T1.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=SUVA, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  ggtitle("Day-10") +
  ylab(expression('SUVA'["254"]~'(mgC'~ L^-1~m^-1*')'))+
  coord_cartesian(ylim = c(0,12))+
  xlab("") +
  Fig.formatting

########## T2
m1.suva.T2.gam=gam(SUVA ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")

m2.suva.T2.gam=gam(SUVA ~  Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")

m3.suva.T2.gam=gam(SUVA ~  s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")

T2.suva.AIC<-AIC(m1.suva.T2.gam, m2.suva.T2.gam, m3.suva.T2.gam)
T2.suva.AIC
##                      df      AIC
## m1.suva.T2.gam 6.560914 166.8072
## m2.suva.T2.gam 5.095612 168.0063
## m3.suva.T2.gam 4.087052 167.3995
summary(m1.suva.T2.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## SUVA ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         4.9241     0.8879   5.546 9.24e-06 ***
## Treatmentunburned  -1.4665     1.2557  -1.168    0.254    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value  
## s(plant.mass..g):Treatmentburned   2.075  2.561 2.763   0.065 .
## s(plant.mass..g):Treatmentunburned 1.000  1.000 0.068   0.796  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.17   Deviance explained = 28.7%
## -REML = 75.261  Scale est. = 11.826    n = 30
T2.suva.anova=anova(m1.suva.T2.gam)
gam.check(m1.suva.T2.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 10 iterations.
## Gradient range [-3.815864e-05,5.391568e-06]
## (score 75.26059 & scale 11.82582).
## Hessian positive definite, eigenvalue range [3.815639e-05,13.02274].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 2.07    1.11    0.65
## s(plant.mass..g):Treatmentunburned 9.00 1.00    1.11    0.56
draw(m1.suva.T2.gam)

concrvity(m1.suva.T2.gam)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5   e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     4.60e-25
## 3 worst    s(plant.mass..g):Treatmentunburned   4.58e-25
## 4 observed para                                 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned     1.54e-30
## 6 observed s(plant.mass..g):Treatmentunburned   9.84e-31
## 7 estimate para                                 5   e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned   1.24e-27
par(mfrow = c(2, 2))
plot(m1.suva.T2.gam, all.terms = TRUE, page=1)

# model predictions
suva.diff.T2<-plot_difference(
  m1.suva.T2.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-5,5))+
  Fig.formatting


###########  
#plot for the model output on rawdata
suva.T2.mod.plot<-
  plot_smooths(
  model = m1.suva.T2.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=SUVA, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  ylim(0, 12) +
  ggtitle("Day-31") +
  ylab("") +
  xlab("plant material (g)") +
  Fig.formatting

suva.T2.mod.plot

########## T3
m1.suva.T3.gam=gam(SUVA ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")

m2.suva.T3.gam=gam(SUVA ~  Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")

m3.suva.T3.gam=gam(SUVA ~  s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")

T3.suva.AIC<-AIC(m1.suva.T3.gam, m2.suva.T3.gam, m3.suva.T3.gam)
T3.suva.AIC
##                       df      AIC
## m1.suva.T3.gam 10.018360 23.87907
## m2.suva.T3.gam  7.067145 28.31202
## m3.suva.T3.gam  5.518278 36.92467
summary(m1.suva.T3.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## SUVA ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        2.95603    0.07731  38.238  < 2e-16 ***
## Treatmentunburned  0.39851    0.10933   3.645  0.00141 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value    
## s(plant.mass..g):Treatmentburned   2.458  3.023 76.33  <2e-16 ***
## s(plant.mass..g):Treatmentunburned 3.270  3.996 92.84  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.955   Deviance explained = 96.5%
## -REML = 14.721  Scale est. = 0.089645  n = 30
T3.suva.anova=anova(m1.suva.T3.gam)
gam.check(m1.suva.T3.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [3.454876e-10,1.6981e-08]
## (score 14.72117 & scale 0.0896449).
## Hessian positive definite, eigenvalue range [0.2269231,13.14457].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 2.46    0.81    0.14
## s(plant.mass..g):Treatmentunburned 9.00 3.27    0.81    0.11
draw(m1.suva.T3.gam)

concrvity(m1.suva.T3.gam)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5   e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     4.60e-25
## 3 worst    s(plant.mass..g):Treatmentunburned   4.58e-25
## 4 observed para                                 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned     4.91e-30
## 6 observed s(plant.mass..g):Treatmentunburned   1.39e-30
## 7 estimate para                                 5   e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned   1.24e-27
par(mfrow = c(2, 2))
plot(m1.suva.T3.gam, all.terms = TRUE, page=1)

# model predictions
suva.diff.T3<-plot_difference(
  m1.suva.T3.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-5,5))+
  Fig.formatting


###########  
#plot for the model output on rawdata
suva.T3.mod.plot<-
  plot_smooths(
  model = m1.suva.T3.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),], 
             aes(x=plant.mass..g, y=SUVA, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  coord_cartesian(ylim = c(0,12))+
  ggtitle("Day-59") +
  ylab("") +
  xlab("") +
  Fig.formatting

########## T4
m1.suva.T4.gam=gam(SUVA ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")

m2.suva.T4.gam=gam(SUVA ~  Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")

m3.suva.T4.gam=gam(SUVA ~  s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")

T4.suva.AIC<-AIC(m1.suva.T4.gam, m2.suva.T4.gam, m3.suva.T4.gam)
T4.suva.AIC
##                      df      AIC
## m1.suva.T4.gam 6.949281 33.36689
## m2.suva.T4.gam 5.560335 38.15182
## m3.suva.T4.gam 4.576257 36.69035
# best is global
summary(m1.suva.T4.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## SUVA ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        2.66735    0.09569  27.874   <2e-16 ***
## Treatmentunburned -0.10136    0.13533  -0.749    0.461    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df      F p-value    
## s(plant.mass..g):Treatmentburned   1.000   1.00 130.62  <2e-16 ***
## s(plant.mass..g):Treatmentunburned 2.531   3.11  52.06  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.909   Deviance explained = 92.3%
## -REML = 17.897  Scale est. = 0.13736   n = 30
T4.suva.anova=anova(m1.suva.T4.gam)
gam.check(m1.suva.T4.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-8.722995e-06,1.789474e-06]
## (score 17.89694 & scale 0.1373614).
## Hessian positive definite, eigenvalue range [8.722846e-06,13.04785].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 1.00    0.98    0.35
## s(plant.mass..g):Treatmentunburned 9.00 2.53    0.98    0.34
draw(m1.suva.T4.gam)

concrvity(m1.suva.T4.gam)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5   e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     4.60e-25
## 3 worst    s(plant.mass..g):Treatmentunburned   4.58e-25
## 4 observed para                                 5   e- 1
## 5 observed s(plant.mass..g):Treatmentburned     5.91e-30
## 6 observed s(plant.mass..g):Treatmentunburned   8.97e-31
## 7 estimate para                                 5   e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned   1.24e-27
par(mfrow = c(2, 2))
plot(m1.suva.T4.gam, all.terms = TRUE, page=1)

# model predictions
suva.diff.T4<-plot_difference(
  m1.suva.T4.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-5,5))+
  Fig.formatting


###########  
#plot for the model output on rawdata
suva.T4.mod.plot<-
  plot_smooths(
  model = m1.suva.T4.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),], 
             aes(x=plant.mass..g, y=SUVA, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  coord_cartesian(ylim = c(0,12))+
  ggtitle("Day-89") +
  ylab("") +
  xlab("") +
  Fig.formatting

suva.T4.mod.plot

mod.suva.rep<-rep(c(
              "Day 10: Treatment with by-factor smooth", 
              "Day 10: Treatment with global smooth", 
              "Day 10: global smooth",
              "Day 31: Treatment with by-factor smooth", 
              "Day 31: Treatment with global smooth", 
              "Day 31: global smooth",
              "Day 59: Treatment with by-factor smooth", 
              "Day 59: Treatment with global smooth", 
              "Day 59: global smooth",
              "Day 89: Treatment with by-factor smooth", 
              "Day 89: Treatment with global smooth", 
              "Day 89: global smooth"))

mod.suva.column<- data.frame(mod.suva.rep)
AIC.suva<-bind_rows(T1.suva.AIC, T2.suva.AIC, T3.suva.AIC, T4.suva.AIC)
AIC.suva.mod<-as.data.frame(cbind(mod.suva.column, AIC.suva))

names(AIC.suva.mod)[1]="Model"


AIC.suva.mod <- AIC.suva.mod %>%
  # Creating an empty column:
  add_column(Variable= c("SUVA"), .before="Model")

write.csv(AIC.suva.mod, "output/AIC models/AIC.suva.csv")

#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/suva.mod.AIC.pdf", width = 6, height = 3.8)       # Export PDF
grid.table(AIC.suva.mod, rows= NULL)
dev.off()
## quartz_off_screen 
##                 2

Generate dataframe for model outputs

####Parametric output
suva.para.sums=as.data.frame(rbind(T1.suva.anova$pTerms.table, T2.suva.anova$pTerms.table, T3.suva.anova$pTerms.table, T4.suva.anova$pTerms.table))
names(suva.para.sums)[1]="df/edf"
suva.para.sums
#create empty column for ref.df
suva.para.sums <- suva.para.sums %>%
  # Creating an empty column:
  add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
  add_column(Ref.df = "-", .after="df/edf")

#####Smooth output
suva.smooth.sums=as.data.frame(rbind(T1.suva.anova$s.table, T2.suva.anova$s.table, T3.suva.anova$s.table, T4.suva.anova$s.table))
names(suva.smooth.sums)[1]="df/edf"
#create empty column for ref.df
suva.smooth.sums <- suva.smooth.sums %>%
  # Creating an empty column:
  add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")

#merge the dataframes
suva.mod.allsums=rbind(suva.para.sums,suva.smooth.sums)
rownames(suva.mod.allsums)<-c(1:12)

#create new ID column
suva.mod.allsums <- suva.mod.allsums %>%
  # Creating an empty column:
  add_column(Variable= c("SUVA"), .before="Model")

#convert dataframe into a table for export
pdf("output/Codys.tables/suva.mod.allsums.pdf", width = 8.5, height = 3.8)       # Export PDF
grid.table(suva.mod.allsums, rows= NULL)
dev.off()

compile raw plots and model-diff plots for final figures


extract.legend <- get_legend(
  # create some space to the left of the legend
  DOC.T1.mod.plot + theme(legend.box.margin = margin(0, 0, 0, 10)))


suva.alltimes<-plot_grid(
  suva.T1.mod.plot+ theme(legend.position = "none"),
  suva.T2.mod.plot+ theme(legend.position = "none"),
  suva.T3.mod.plot+ theme(legend.position = "none"),
  suva.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
  rel_widths = c(8,8,8,8,3), ncol=5)

suva.alltimes
ggsave("figures/SUVA.alltimes.mod.pdf", height=4, width=12)

suva.mod.diffs<-plot_grid(
  suva.diff.T1+ theme(legend.position = "none")+ ggtitle("SUVA - Day-10"),
  suva.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
  suva.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
  suva.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
  rel_widths = c(8,8,8,8), ncol=4)

suva.mod.diffs
ggsave("figures/suva.mod.diffs.pdf", height=4, width=12)

Sr analysis and plots

######## T0 model
m1.sr.T0.gam=gam(sr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")

m2.sr.T0.gam=gam(sr ~  Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")

m3.sr.T0.gam=gam(sr ~  s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")

T0.sr.AIC<-AIC(m1.sr.T0.gam, m2.sr.T0.gam, m3.sr.T0.gam)
T0.sr.AIC
##                    df      AIC
## m1.sr.T0.gam 5.656040 72.39828
## m2.sr.T0.gam 4.000309 76.34838
## m3.sr.T0.gam 3.000189 77.53584
summary(m1.sr.T0.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## sr ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         3.8530     0.1871  20.591   <2e-16 ***
## Treatmentunburned  -0.5062     0.2646  -1.913    0.067 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value  
## s(plant.mass..g):Treatmentburned   1.000  1.000 6.896  0.0143 *
## s(plant.mass..g):Treatmentunburned 1.376  1.656 0.401  0.5296  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.234   Deviance explained = 32.3%
## -REML = 34.155  Scale est. = 0.52519   n = 30
anova.gam(m1.sr.T0.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## sr ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 3.659   0.067
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   1.000  1.000 6.896  0.0143
## s(plant.mass..g):Treatmentunburned 1.376  1.656 0.401  0.5296
gam.check(m1.sr.T0.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-1.375453e-05,2.556971e-06]
## (score 34.15521 & scale 0.5251928).
## Hessian positive definite, eigenvalue range [1.375414e-05,13.00273].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 1.00     0.9    0.19
## s(plant.mass..g):Treatmentunburned 9.00 1.38     0.9    0.24
draw(m1.sr.T0.gam)

concrvity(m1.sr.T0.gam)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5   e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     4.60e-25
## 3 worst    s(plant.mass..g):Treatmentunburned   4.58e-25
## 4 observed para                                 5   e- 1
## 5 observed s(plant.mass..g):Treatmentburned     5.91e-30
## 6 observed s(plant.mass..g):Treatmentunburned   9.79e-31
## 7 estimate para                                 5   e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned   1.24e-27
par(mfrow = c(2, 2))
plot(m1.sr.T0.gam, all.terms = TRUE, page=1)

# model predictions
sr.diff.T0<-plot_difference(
  m1.sr.T0.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-0.75,0.75))+
  Fig.formatting

###########  
sr.T0.mod.plot<-
  plot_smooths(
  model = m3.sr.T0.gam,
  series = plant.mass..g,
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),], 
             aes(x=plant.mass..g, y=sr, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  ggtitle("Day-0") +
  ylab(expression('S'[italic(R)]))+
  xlab("plant material (g)") +
  Fig.formatting

sr.T0.mod.plot

######## T1 model
m1.sr.T1.gam=gam(sr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")

m2.sr.T1.gam=gam(sr ~  Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")

m3.sr.T1.gam=gam(sr ~  s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")

T1.sr.AIC<-AIC(m1.sr.T1.gam, m2.sr.T1.gam, m3.sr.T1.gam)
T1.sr.AIC
##                    df       AIC
## m1.sr.T1.gam 7.127378  7.775483
## m2.sr.T1.gam 5.172277 12.562474
## m3.sr.T1.gam 4.186417 11.180366
summary(m1.sr.T1.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## sr ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        1.29977    0.06211   20.93   <2e-16 ***
## Treatmentunburned  0.07200    0.08784    0.82     0.42    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F  p-value    
## s(plant.mass..g):Treatmentburned   1.000  1.000 12.85  0.00149 ** 
## s(plant.mass..g):Treatmentunburned 2.545  3.127 13.36 2.25e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.643   Deviance explained = 69.9%
## -REML = 6.6791  Scale est. = 0.057874  n = 30
T1.sr.anova=anova(m1.sr.T1.gam)
gam.check(m1.sr.T1.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 12 iterations.
## Gradient range [-2.139521e-06,2.098115e-07]
## (score 6.679058 & scale 0.05787405).
## Hessian positive definite, eigenvalue range [2.139513e-06,13.04761].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value   
## s(plant.mass..g):Treatmentburned   9.00 1.00    0.65   0.015 * 
## s(plant.mass..g):Treatmentunburned 9.00 2.55    0.65   0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
draw(m1.sr.T1.gam)

concrvity(m1.sr.T1.gam)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5   e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     4.60e-25
## 3 worst    s(plant.mass..g):Treatmentunburned   4.58e-25
## 4 observed para                                 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned     5.91e-30
## 6 observed s(plant.mass..g):Treatmentunburned   7.24e-31
## 7 estimate para                                 5   e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned   1.24e-27
par(mfrow = c(2, 2))
plot(m1.sr.T1.gam, all.terms = TRUE, page=1)

# model predictions
sr.diff.T1<-plot_difference(
  m1.sr.T1.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-0.75,0.75))+
  Fig.formatting

###########  
#plot for the model output on rawdata
sr.T1.mod.plot<-
  plot_smooths(
  model = m1.sr.T1.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=sr, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  ggtitle("Day-10") +
  ylab(expression('S'[italic(R)]))+
  coord_cartesian(ylim = c(0.25,2.5))+
  xlab("") +
  Fig.formatting

sr.T1.mod.plot

########## T2
m1.sr.T2.gam=gam(sr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")

m2.sr.T2.gam=gam(sr ~  Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")

m3.sr.T2.gam=gam(sr ~  s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")

T2.sr.AIC<-AIC(m1.sr.T2.gam, m2.sr.T2.gam, m3.sr.T2.gam)
T2.sr.AIC
##                     df       AIC
## m1.sr.T2.gam 11.537024 -13.51026
## m2.sr.T2.gam  7.677441 -20.11997
## m3.sr.T2.gam  6.712723 -22.00780
summary(m3.sr.T2.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## sr ~ s(plant.mass..g)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.4040     0.0269   52.19   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df     F  p-value    
## s(plant.mass..g) 4.171  5.053 17.25 8.53e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.748   Deviance explained = 78.4%
## -REML = -6.5014  Scale est. = 0.021714  n = 30
T2.sr.anova=anova(m1.sr.T2.gam)
gam.check(m3.sr.T2.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-2.027443e-09,5.61081e-10]
## (score -6.501377 & scale 0.02171413).
## Hessian positive definite, eigenvalue range [1.070181,14.19444].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                    k'  edf k-index p-value
## s(plant.mass..g) 9.00 4.17    1.59       1
draw(m3.sr.T2.gam)

concrvity(m3.sr.T2.gam)
## # A tibble: 6 × 3
##   type     term             concurvity
##   <chr>    <chr>                 <dbl>
## 1 worst    para               4.61e-25
## 2 worst    s(plant.mass..g)   4.61e-25
## 3 observed para               4.61e-25
## 4 observed s(plant.mass..g)   3.94e-31
## 5 estimate para               4.61e-25
## 6 estimate s(plant.mass..g)   1.24e-27
par(mfrow = c(2, 2))
plot(m3.sr.T2.gam, all.terms = TRUE, page=1)

# model predictions
sr.diff.T2<-plot_difference(
  m1.sr.T2.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-0.75,0.75))+
  Fig.formatting

###########  
#plot for the model output on rawdata
sr.T2.mod.plot<-
  plot_smooths(
  model = m3.sr.T2.gam,
  series = plant.mass..g,
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=sr, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  coord_cartesian(ylim = c(0.25,2.5))+
  ggtitle("Day-31") +
  ylab("") +
  xlab("plant material (g)") +
  Fig.formatting

sr.T2.mod.plot 

########## T3
m1.sr.T3.gam=gam(sr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")

m2.sr.T3.gam=gam(sr ~  Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")

m3.sr.T3.gam=gam(sr ~  s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")

m4.sr.T3.gam=gam(sr ~  Treatment + s(plant.mass..g, m=c(2,0)), subset = Time.point=="T3", data = DOC.df, method = "ML")

T3.sr.AIC<-AIC(m1.sr.T3.gam, m2.sr.T3.gam, m3.sr.T3.gam)
T3.sr.AIC
##                    df       AIC
## m1.sr.T3.gam 9.688689 -1.124062
## m2.sr.T3.gam 6.776426 -6.818488
## m3.sr.T3.gam 5.747507 -6.468986
summary(m1.sr.T3.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## sr ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        1.39164    0.05136  27.098   <2e-16 ***
## Treatmentunburned -0.09869    0.07263  -1.359    0.188    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F  p-value    
## s(plant.mass..g):Treatmentburned   3.035  3.715 11.99 4.06e-05 ***
## s(plant.mass..g):Treatmentunburned 2.548  3.131 11.13 0.000108 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.729   Deviance explained =   79%
## -REML = 3.8444  Scale est. = 0.039562  n = 30
T3.sr.anova=anova(m1.sr.T3.gam)
gam.check(m1.sr.T3.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-1.513811e-11,1.06315e-12]
## (score 3.844409 & scale 0.03956209).
## Hessian positive definite, eigenvalue range [0.6287135,13.13229].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 3.04    1.05    0.52
## s(plant.mass..g):Treatmentunburned 9.00 2.55    1.05    0.50
draw(m1.sr.T3.gam)

concrvity(m1.sr.T3.gam)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5   e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     4.60e-25
## 3 worst    s(plant.mass..g):Treatmentunburned   4.58e-25
## 4 observed para                                 5   e- 1
## 5 observed s(plant.mass..g):Treatmentburned     1.64e-30
## 6 observed s(plant.mass..g):Treatmentunburned   6.91e-31
## 7 estimate para                                 5   e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned   1.24e-27
par(mfrow = c(2, 2))
plot(m1.sr.T3.gam, all.terms = TRUE, page=1)

# model predictions
sr.diff.T3<-plot_difference(
  m1.sr.T3.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-0.75,0.75))+
  Fig.formatting

###########  
#plot for the model output on rawdata
sr.T3.mod.plot<-
  plot_smooths(
  model = m2.sr.T3.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),], 
             aes(x=plant.mass..g, y=sr, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  scale_fill_manual(values = c("brown1", "mediumseagreen"))+
  coord_cartesian(ylim = c(0.25,2.5))+
  geom_line(aes(linetype=Treatment))+
  ggtitle("Day-59") +
  ylab("") +
  xlab("") +
  Fig.formatting

sr.T3.mod.plot

########## T4
m1.sr.T4.gam=gam(sr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")

m2.sr.T4.gam=gam(sr ~  Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")

m3.sr.T4.gam=gam(sr ~  s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")

T4.sr.AIC<-AIC(m1.sr.T4.gam, m2.sr.T4.gam, m3.sr.T4.gam)
T4.sr.AIC
##                     df       AIC
## m1.sr.T4.gam 11.398371 -57.07811
## m2.sr.T4.gam  7.363534 -42.13080
## m3.sr.T4.gam  6.396737 -44.12181
summary(m1.sr.T4.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## sr ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.288325   0.019910  64.709   <2e-16 ***
## Treatmentunburned 0.004562   0.028156   0.162    0.873    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value    
## s(plant.mass..g):Treatmentburned   2.890  3.542 25.64  <2e-16 ***
## s(plant.mass..g):Treatmentunburned 4.498  5.427 46.63  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.922   Deviance explained = 94.4%
## -REML = -17.839  Scale est. = 0.0059459  n = 30
T4.sr.anova=anova(m1.sr.T4.gam)
gam.check(m1.sr.T4.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [-5.197224e-09,4.086935e-10]
## (score -17.83887 & scale 0.005945903).
## Hessian positive definite, eigenvalue range [0.6347134,13.32902].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 2.89       1    0.45
## s(plant.mass..g):Treatmentunburned 9.00 4.50       1    0.47
draw(m1.sr.T4.gam)

concrvity(m1.sr.T4.gam)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5   e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     4.60e-25
## 3 worst    s(plant.mass..g):Treatmentunburned   4.58e-25
## 4 observed para                                 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned     3.17e-30
## 6 observed s(plant.mass..g):Treatmentunburned   9.19e-31
## 7 estimate para                                 5   e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned   1.24e-27
par(mfrow = c(2, 2))
plot(m1.sr.T4.gam, all.terms = TRUE, page=1)

# model predictions
sr.diff.T4<-plot_difference(
  m1.sr.T4.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-0.75,0.75))+
  Fig.formatting

###########  
#plot for the model output on rawdata
sr.T4.mod.plot<-
  plot_smooths(
  model = m1.sr.T4.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),], 
             aes(x=plant.mass..g, y=sr, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  coord_cartesian(ylim = c(0.25,2.5))+
  ggtitle("Day-89") +
  ylab("") +
  xlab("") +
  Fig.formatting

sr.T4.mod.plot

mod.sr.rep<-rep(c(
              "Day 10: Treatment with by-factor smooth", 
              "Day 10: Treatment with global smooth", 
              "Day 10: global smooth",
              "Day 31: Treatment with by-factor smooth", 
              "Day 31: Treatment with global smooth", 
              "Day 31: global smooth",
              "Day 59: Treatment with by-factor smooth", 
              "Day 59: Treatment with global smooth", 
              "Day 59: global smooth",
              "Day 89: Treatment with by-factor smooth", 
              "Day 89: Treatment with global smooth", 
              "Day 89: global smooth"))

mod.sr.column<- data.frame(mod.sr.rep)
AIC.sr<-bind_rows(T1.sr.AIC, T2.sr.AIC, T3.sr.AIC, T4.sr.AIC)
AIC.sr
##                     df        AIC
## m1.sr.T1.gam  7.127378   7.775483
## m2.sr.T1.gam  5.172277  12.562474
## m3.sr.T1.gam  4.186417  11.180366
## m1.sr.T2.gam 11.537024 -13.510260
## m2.sr.T2.gam  7.677441 -20.119965
## m3.sr.T2.gam  6.712723 -22.007797
## m1.sr.T3.gam  9.688689  -1.124062
## m2.sr.T3.gam  6.776426  -6.818488
## m3.sr.T3.gam  5.747507  -6.468986
## m1.sr.T4.gam 11.398371 -57.078114
## m2.sr.T4.gam  7.363534 -42.130799
## m3.sr.T4.gam  6.396737 -44.121805
AIC.sr.mod<-as.data.frame(cbind(mod.sr.column, AIC.sr))

names(AIC.sr.mod)[1]="Model"

AIC.sr.mod <- AIC.sr.mod %>%
  # Creating an empty column:
  add_column(Variable= c("Slope ratio"), .before="Model")

write.csv(AIC.sr.mod, "output/AIC models/AIC.sr.csv")

#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/sr.mod.AIC.pdf", width = 6.3, height = 3.8)       # Export PDF
grid.table(AIC.sr.mod, rows=NULL)
dev.off()
## quartz_off_screen 
##                 2

Generate dataframe for model outputs

####Parametric output
sr.para.sums=as.data.frame(rbind(T1.sr.anova$pTerms.table, T2.sr.anova$pTerms.table, T3.sr.anova$pTerms.table, T4.sr.anova$pTerms.table))
names(sr.para.sums)[1]="df/edf"
#create empty column for ref.df
sr.para.sums <- sr.para.sums %>%
  # Creating an empty column:
  add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
  add_column(Ref.df = "-", .after="df/edf")
  

#####Smooth output
sr.smooth.sums=as.data.frame(rbind(T1.sr.anova$s.table, T2.sr.anova$s.table, T3.sr.anova$s.table, T4.sr.anova$s.table))
names(sr.smooth.sums)[1]="df/edf"
#create empty column for ref.df
sr.smooth.sums <- sr.smooth.sums %>%
  # Creating an empty column:
  add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")

sr.smooth.sums
#merge the dataframes
sr.mod.allsums=rbind(sr.para.sums,sr.smooth.sums)
rownames(sr.mod.allsums)<-c(1:12)

#create new ID column
sr.mod.allsums <- sr.mod.allsums %>%
  # Creating an empty column:
  add_column(Variable= c("Slope ratio"), .before="Model")

#convert dataframe into a table for export
pdf("output/Codys.tables/sr.mod.allsums.pdf", width = 8.5, height = 3.8)       # Export PDF
grid.table(sr.mod.allsums, rows= NULL)
dev.off()

Sr: compile raw plots and model-diff plots for final figures

sr.alltimes<-plot_grid(
  sr.T1.mod.plot+ theme(legend.position = "none"),
  sr.T2.mod.plot+ theme(legend.position = "none"),
  sr.T3.mod.plot+ theme(legend.position = "none"),
  sr.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
  rel_widths = c(8,8,8,8,3), ncol=5)

sr.alltimes

ggsave("figures/sr.alltimes.mod.pdf", height=4, width=12)

sr.mod.diffs<-plot_grid(
  sr.diff.T1+ theme(legend.position = "none")+ ggtitle("Sr - Day-10"),
  sr.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
  sr.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
  sr.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
  rel_widths = c(8,8,8,8), ncol=4)

sr.mod.diffs

ggsave("figures/sr.mod.diffs.pdf", height=4, width=12)

HIX analysis and plots

######## T0 model
m1.hix.T0.gam=gam(hix ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")

m2.hix.T0.gam=gam(hix ~  Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")

m3.hix.T0.gam=gam(hix ~  s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")

T0.hix.AIC<-AIC(m1.hix.T0.gam, m2.hix.T0.gam, m3.hix.T0.gam)
T0.hix.AIC
##                     df      AIC
## m1.hix.T0.gam 5.700361 20.37557
## m2.hix.T0.gam 5.240725 18.48445
## m3.hix.T0.gam 3.985728 23.07642
summary(m1.hix.T0.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## hix ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        1.52625    0.07855  19.429   <2e-16 ***
## Treatmentunburned -0.27629    0.11109  -2.487   0.0197 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   1.406    1.7 0.489   0.693
## s(plant.mass..g):Treatmentunburned 1.000    1.0 2.258   0.145
## 
## R-sq.(adj) =  0.178   Deviance explained = 27.4%
## -REML = 11.609  Scale est. = 0.092562  n = 30
anova.gam(m1.hix.T0.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## hix ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 6.185  0.0197
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned   1.406  1.700 0.489   0.693
## s(plant.mass..g):Treatmentunburned 1.000  1.000 2.258   0.145
gam.check(m1.hix.T0.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 9 iterations.
## Gradient range [-6.467201e-06,1.475016e-05]
## (score 11.60873 & scale 0.09256152).
## Hessian positive definite, eigenvalue range [6.46737e-06,13.00317].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 1.41    0.95    0.39
## s(plant.mass..g):Treatmentunburned 9.00 1.00    0.95    0.28
draw(m1.hix.T0.gam)

concrvity(m1.hix.T0.gam)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5   e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     4.60e-25
## 3 worst    s(plant.mass..g):Treatmentunburned   4.58e-25
## 4 observed para                                 5   e- 1
## 5 observed s(plant.mass..g):Treatmentburned     9.13e-30
## 6 observed s(plant.mass..g):Treatmentunburned   9.84e-31
## 7 estimate para                                 5   e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned   1.24e-27
par(mfrow = c(2, 2))
plot(m1.hix.T0.gam, all.terms = TRUE, page=1)

# model predictions
hix.diff.T0<-plot_difference(
  m1.hix.T0.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
hix.T0.mod.plot<-
  plot_smooths(
  model = m3.hix.T0.gam,
  series = plant.mass..g,
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),], 
             aes(x=plant.mass..g, y=hix, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  ggtitle("Day-0") +
    coord_cartesian(ylim = c(0,6))+
  ylab("HIX")+
  xlab("plant material (g)") +
  Fig.formatting


######## T1 model
m1.hix.T1.gam=gam(hix ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")

m2.hix.T1.gam=gam(hix ~  Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")

m3.hix.T1.gam=gam(hix ~  s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")

T1.hix.AIC<-AIC(m1.hix.T1.gam, m2.hix.T1.gam, m3.hix.T1.gam)
T1.hix.AIC
##                     df      AIC
## m1.hix.T1.gam 5.000061 62.10694
## m2.hix.T1.gam 4.000060 60.23250
## m3.hix.T1.gam 3.000061 59.42114
summary(m3.hix.T1.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## hix ~ s(plant.mass..g)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.7544     0.1114   15.75 1.92e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                  edf Ref.df     F  p-value    
## s(plant.mass..g)   1      1 15.64 0.000475 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.335   Deviance explained = 35.8%
## -REML = 29.297  Scale est. = 0.37225   n = 30
T1.hix.anova=anova(m1.hix.T1.gam)
gam.check(m3.sr.T1.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-6.202238e-07,-3.062326e-07]
## (score 6.592856 & scale 0.07082567).
## Hessian positive definite, eigenvalue range [0.2360873,14.01068].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                    k'  edf k-index p-value  
## s(plant.mass..g) 9.00 1.77    0.73   0.035 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
draw(m3.hix.T1.gam)

concrvity(m3.hix.T1.gam)
## # A tibble: 6 × 3
##   type     term             concurvity
##   <chr>    <chr>                 <dbl>
## 1 worst    para               4.61e-25
## 2 worst    s(plant.mass..g)   4.61e-25
## 3 observed para               4.61e-25
## 4 observed s(plant.mass..g)   2.60e-32
## 5 estimate para               4.61e-25
## 6 estimate s(plant.mass..g)   1.24e-27
par(mfrow = c(2, 2))
plot(m3.hix.T1.gam, all.terms = TRUE, page=1)

# model predictions
hix.diff.T1<-plot_difference(
  m1.hix.T1.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-2.5,1.5))+
  Fig.formatting

###########  
#plot for the model output on rawdata
hix.T1.mod.plot<-
  plot_smooths(
  model = m3.hix.T1.gam,
  series = plant.mass..g,
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=hix, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  coord_cartesian(ylim = c(0,6))+
  ggtitle("Day-10") +
  ylab("HIX")+
  xlab("") +
  Fig.formatting
hix.T1.mod.plot

########## T2
m1.hix.T2.gam=gam(hix ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")

m2.hix.T2.gam=gam(hix ~  Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")

m3.hix.T2.gam=gam(hix ~  s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")

m4.hix.T2.gam=gam(hix ~  Treatment + s(plant.mass..g, m=c(2,0)), subset = Time.point=="T2", data = DOC.df, method = "REML")

T2.hix.AIC<-AIC(m1.hix.T2.gam, m2.hix.T2.gam, m3.hix.T2.gam)
T2.hix.AIC
##                     df      AIC
## m1.hix.T2.gam 5.000022 30.35861
## m2.hix.T2.gam 4.000030 28.41308
## m3.hix.T2.gam 3.000045 28.88814
summary(m2.hix.T2.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## hix ~ Treatment + s(plant.mass..g)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        2.17247    0.09254  23.475   <2e-16 ***
## Treatmentunburned  0.19943    0.13088   1.524    0.139    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                  edf Ref.df     F p-value  
## s(plant.mass..g)   1      1 6.651  0.0157 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.194   Deviance explained = 24.9%
## -REML = 15.017  Scale est. = 0.12847   n = 30
T2.hix.anova=anova(m1.hix.T2.gam)
gam.check(m2.hix.T2.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 13 iterations.
## Gradient range [-5.861103e-06,1.733725e-06]
## (score 15.01688 & scale 0.1284671).
## Hessian positive definite, eigenvalue range [5.861019e-06,13.5].
## Model rank =  11 / 11 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                  k' edf k-index p-value
## s(plant.mass..g)  9   1    0.86    0.12
draw(m2.hix.T2.gam)

concrvity(m2.hix.T2.gam)
## # A tibble: 6 × 3
##   type     term             concurvity
##   <chr>    <chr>                 <dbl>
## 1 worst    para               5   e- 1
## 2 worst    s(plant.mass..g)   4.61e-25
## 3 observed para               5   e- 1
## 4 observed s(plant.mass..g)   2.96e-32
## 5 estimate para               5   e- 1
## 6 estimate s(plant.mass..g)   1.24e-27
par(mfrow = c(2, 2))
plot(m2.hix.T2.gam, all.terms = TRUE, page=1)

# model predictions
hix.diff.T2<-plot_difference(
  m1.hix.T2.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-2.5,1.5))+
  Fig.formatting

###########  
#plot for the model output on rawdata
hix.T2.mod.plot<-
  plot_smooths(
  model = m2.hix.T2.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=hix, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  coord_cartesian(ylim = c(0,6))+
  geom_line(aes(linetype=Treatment))+
  ggtitle("Day-31") +
  ylab("") +
  xlab("plant material (g)") +
  Fig.formatting

hix.T2.mod.plot

########## T3
m1.hix.T3.gam=gam(hix ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")

m2.hix.T3.gam=gam(hix ~  Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")

m3.hix.T3.gam=gam(hix ~  s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")

T3.hix.AIC<-AIC(m1.hix.T3.gam, m2.hix.T3.gam, m3.hix.T3.gam)
T3.hix.AIC
##                      df      AIC
## m1.hix.T3.gam 11.953543 23.11911
## m2.hix.T3.gam  7.757426 23.50220
## m3.hix.T3.gam  5.936863 35.52011
summary(m1.hix.T3.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## hix ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         2.4393     0.0745  32.741  < 2e-16 ***
## Treatmentunburned   0.4249     0.1054   4.033 0.000623 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F  p-value    
## s(plant.mass..g):Treatmentburned   2.841  3.484 7.129 0.001285 ** 
## s(plant.mass..g):Treatmentunburned 4.609  5.553 6.346 0.000791 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.718   Deviance explained =   80%
## -REML = 16.611  Scale est. = 0.083257  n = 30
T3.hix.anova=anova(m1.hix.T3.gam)
gam.check(m1.hix.T3.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-1.41462e-08,1.860842e-09]
## (score 16.61062 & scale 0.08325724).
## Hessian positive definite, eigenvalue range [0.6660674,13.3327].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 2.84    0.91    0.26
## s(plant.mass..g):Treatmentunburned 9.00 4.61    0.91    0.23
draw(m1.hix.T3.gam)

concrvity(m1.hix.T3.gam)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5   e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     4.60e-25
## 3 worst    s(plant.mass..g):Treatmentunburned   4.58e-25
## 4 observed para                                 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned     9.86e-31
## 6 observed s(plant.mass..g):Treatmentunburned   2.60e-30
## 7 estimate para                                 5   e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned   1.24e-27
par(mfrow = c(2, 2))
plot(m1.hix.T3.gam, all.terms = TRUE, page=1)

# model predictions
hix.diff.T3<-plot_difference(
  m1.hix.T3.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-2.5,1.5))+
  Fig.formatting

###########  
#plot for the model output on rawdata
hix.T3.mod.plot<-
  plot_smooths(
  model = m1.hix.T3.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),], 
             aes(x=plant.mass..g, y=hix, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  coord_cartesian(ylim = c(0,6))+
  ggtitle("Day-59") +
  ylab("") +
  xlab("") +
  Fig.formatting

########## T4
m1.hix.T4.gam=gam(hix ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")

m2.hix.T4.gam=gam(hix ~  Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")

m3.hix.T4.gam=gam(hix ~  s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")

T4.hix.AIC<-AIC(m1.hix.T4.gam, m2.hix.T4.gam, m3.hix.T4.gam)
T4.hix.AIC
##                      df      AIC
## m1.hix.T4.gam 10.053318 25.12183
## m2.hix.T4.gam  6.208475 41.21670
## m3.hix.T4.gam  4.841730 52.43251
summary(m1.hix.T4.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## hix ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        2.59052    0.07888  32.840  < 2e-16 ***
## Treatmentunburned  0.58969    0.11156   5.286 2.55e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df     F  p-value    
## s(plant.mass..g):Treatmentburned   2.501  3.074 10.75 0.000141 ***
## s(plant.mass..g):Treatmentunburned 3.257  3.980 37.18  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.877   Deviance explained = 90.5%
## -REML =  15.28  Scale est. = 0.09334   n = 30
T4.hix.anova=anova(m1.hix.T4.gam)
gam.check(m1.hix.T4.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 3 iterations.
## Gradient range [-6.675986e-06,2.083557e-06]
## (score 15.28009 & scale 0.0933395).
## Hessian positive definite, eigenvalue range [0.3536202,13.14571].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value  
## s(plant.mass..g):Treatmentburned   9.00 2.50    0.73   0.070 .
## s(plant.mass..g):Treatmentunburned 9.00 3.26    0.73   0.045 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
draw(m1.hix.T4.gam)

concrvity(m1.hix.T4.gam)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5   e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     4.60e-25
## 3 worst    s(plant.mass..g):Treatmentunburned   4.58e-25
## 4 observed para                                 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned     1.95e-30
## 6 observed s(plant.mass..g):Treatmentunburned   2.31e-30
## 7 estimate para                                 5   e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned   1.24e-27
par(mfrow = c(2, 2))
plot(m1.hix.T4.gam, all.terms = TRUE, page=1)

# model predictions
hix.diff.T4<-plot_difference(
  m1.hix.T4.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-2.5,1.5))+
  Fig.formatting

###########  
#plot for the model output on rawdata
hix.T4.mod.plot<-
  plot_smooths(
  model = m1.hix.T4.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),], 
             aes(x=plant.mass..g, y=hix, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  coord_cartesian(ylim = c(0,6))+
  ggtitle("Day-89") +
  ylab("") +
  xlab("") +
  Fig.formatting

mod.hix.rep<-rep(c(
              "Day 10: Treatment with by-factor smooth", 
              "Day 10: Treatment with global smooth", 
              "Day 10: global smooth",
              "Day 31: Treatment with by-factor smooth", 
              "Day 31: Treatment with global smooth", 
              "Day 31: global smooth",
              "Day 59: Treatment with by-factor smooth", 
              "Day 59: Treatment with global smooth", 
              "Day 59: global smooth",
              "Day 89: Treatment with by-factor smooth", 
              "Day 89: Treatment with global smooth", 
              "Day 89: global smooth"))



mod.hix.column<- data.frame(mod.hix.rep)
AIC.hix<-bind_rows(T1.hix.AIC, T2.hix.AIC, T3.hix.AIC, T4.hix.AIC)

AIC.hix.mod<-as.data.frame(cbind(mod.hix.column, AIC.hix))

names(AIC.hix.mod)[1]="Model"

AIC.hix.mod <- AIC.hix.mod %>%
  # Creating an empty column:
  add_column(Variable= c("HIX"), .before="Model")


write.csv(AIC.hix.mod, "output/AIC models/AIC.hix.csv")

#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/hix.mod.AIC.pdf", width = 6, height = 3.8)       # Export PDF
grid.table(AIC.hix.mod, rows= NULL)
dev.off()
## quartz_off_screen 
##                 2

Generate dataframe for model outputs

####Parametric output
hix.para.sums=as.data.frame(rbind(T1.hix.anova$pTerms.table, T2.hix.anova$pTerms.table, T3.hix.anova$pTerms.table, T4.hix.anova$pTerms.table))
names(hix.para.sums)[1]="df/edf"
#create empty column for ref.df
hix.para.sums <- hix.para.sums %>%
  # Creating an empty column:
  add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
  add_column(Ref.df = "-", .after="df/edf")
  

#####Smooth output
hix.smooth.sums=as.data.frame(rbind(T1.hix.anova$s.table, T2.hix.anova$s.table, T3.hix.anova$s.table, T4.hix.anova$s.table))
names(hix.smooth.sums)[1]="df/edf"
#create empty column for ref.df
hix.smooth.sums <- hix.smooth.sums %>%
  # Creating an empty column:
  add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")


#merge the dataframes
hix.mod.allsums=rbind(hix.para.sums,hix.smooth.sums)


#create new ID column
hix.mod.allsums <- hix.mod.allsums %>%
  # Creating an empty column:
  add_column(Variable= c("HIX"), .before="Model")


#convert dataframe into a table for export
pdf("output/Codys.tables/hix.mod.allsums.pdf", width = 8.3, height = 3.8)       # Export PDF
grid.table(hix.mod.allsums, rows= NULL)
dev.off()
## quartz_off_screen 
##                 2

HIX: compile raw plots and model-diff plots for final figures

hix.alltimes<-plot_grid(
  hix.T1.mod.plot+ theme(legend.position = "none"),
    hix.T2.mod.plot+ theme(legend.position = "none"),
   hix.T3.mod.plot+ theme(legend.position = "none"),
    hix.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
  rel_widths = c(8,8,8,8,3), ncol=5)

hix.alltimes

ggsave("figures/hix.alltimes.mod.pdf", height=4, width=12)

hix.mod.diffs<-plot_grid(
  hix.diff.T1+ theme(legend.position = "none")+ ggtitle("HIX - Day-10"),
  hix.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
  hix.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
  hix.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
  rel_widths = c(8,8,8,8), ncol=4)

hix.mod.diffs

ggsave("figures/hix.mod.diffs.pdf", height=4, width=12)

Freshness analysis and plots

######## T0 model
m1.fr.T0.gam=gam(fr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")

m2.fr.T0.gam=gam(fr ~  Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")

m3.fr.T0.gam=gam(fr ~  s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")

T0.fr.AIC<-AIC(m1.fr.T0.gam, m2.fr.T0.gam, m3.fr.T0.gam)
T0.fr.AIC
##                    df       AIC
## m1.fr.T0.gam 5.000114 -56.31473
## m2.fr.T0.gam 4.000035 -56.48287
## m3.fr.T0.gam 3.000036 -58.25380
summary(m1.fr.T0.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## fr ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.89630    0.02222   40.33   <2e-16 ***
## Treatmentunburned  0.01446    0.03143    0.46    0.649    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                    edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned     1      1 1.038   0.318
## s(plant.mass..g):Treatmentunburned   1      1 0.625   0.436
## 
## R-sq.(adj) =  -0.0404   Deviance explained = 6.73%
## -REML = -21.46  Scale est. = 0.0074074  n = 30
anova.gam(m1.fr.T0.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## fr ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric Terms:
##           df     F p-value
## Treatment  1 0.212   0.649
## 
## Approximate significance of smooth terms:
##                                    edf Ref.df     F p-value
## s(plant.mass..g):Treatmentburned     1      1 1.038   0.318
## s(plant.mass..g):Treatmentunburned   1      1 0.625   0.436
gam.check(m1.fr.T0.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 12 iterations.
## Gradient range [-1.342278e-05,1.476944e-05]
## (score -21.46003 & scale 0.007407422).
## Hessian positive definite, eigenvalue range [4.232612e-07,12.99999].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                    k' edf k-index p-value
## s(plant.mass..g):Treatmentburned    9   1    1.13    0.75
## s(plant.mass..g):Treatmentunburned  9   1    1.13    0.72
draw(m1.fr.T0.gam)

concrvity(m1.fr.T0.gam)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5   e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     4.60e-25
## 3 worst    s(plant.mass..g):Treatmentunburned   4.58e-25
## 4 observed para                                 5.00e- 1
## 5 observed s(plant.mass..g):Treatmentburned     5.91e-30
## 6 observed s(plant.mass..g):Treatmentunburned   9.84e-31
## 7 estimate para                                 5   e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned   1.24e-27
par(mfrow = c(2, 2))
plot(m1.fr.T0.gam, all.terms = TRUE, page=1)

# model predictions
fr.diff.T0<-plot_difference(
  m1.fr.T0.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)

###########  
fr.T0.mod.plot<-
  plot_smooths(
  model = m3.fr.T0.gam,
  series = plant.mass..g,
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),], 
             aes(x=plant.mass..g, y=fr, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  ggtitle("Day-0") + 
  ylim(0.25,1.75)+
  xlab("plant material (g)") +
  Fig.formatting

fr.T0.mod.plot

######## T1 model
m1.fr.T1.gam=gam(fr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")

m2.fr.T1.gam=gam(fr ~  Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")

m3.fr.T1.gam=gam(fr ~  s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")

T1.fr.AIC<-AIC(m1.fr.T1.gam, m2.fr.T1.gam, m3.fr.T1.gam)
T1.fr.AIC
##                    df        AIC
## m1.fr.T1.gam 6.194913 -5.4939710
## m2.fr.T1.gam 4.000010 -0.6819142
## m3.fr.T1.gam 3.000009  2.7456950
summary(m1.fr.T1.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## fr ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.67640    0.05057  13.375 5.89e-13 ***
## Treatmentunburned  0.18648    0.07152   2.607   0.0151 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df      F  p-value    
## s(plant.mass..g):Treatmentburned   1.774  2.195  0.971 0.446767    
## s(plant.mass..g):Treatmentunburned 1.000  1.000 18.832 0.000207 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.46   Deviance explained =   53%
## -REML = 0.45251  Scale est. = 0.038362  n = 30
T1.fr.anova=anova(m1.fr.T1.gam)
gam.check(m1.fr.T1.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-8.295506e-07,2.641631e-06]
## (score 0.4525131 & scale 0.0383616).
## Hessian positive definite, eigenvalue range [8.295565e-07,13.01166].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 1.77     0.9    0.22
## s(plant.mass..g):Treatmentunburned 9.00 1.00     0.9    0.14
draw(m1.fr.T1.gam)

concrvity(m1.fr.T1.gam)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5   e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     4.60e-25
## 3 worst    s(plant.mass..g):Treatmentunburned   4.58e-25
## 4 observed para                                 5   e- 1
## 5 observed s(plant.mass..g):Treatmentburned     6.67e-31
## 6 observed s(plant.mass..g):Treatmentunburned   9.84e-31
## 7 estimate para                                 5   e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned   1.24e-27
par(mfrow = c(2, 2))
plot(m1.fr.T1.gam, all.terms = TRUE, page=1)

# model predictions
fr.diff.T1<-plot_difference(
  m1.fr.T1.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-1,0.2))+
  Fig.formatting

###########  
#plot for the model output on rawdata
fr.T1.mod.plot<-
  plot_smooths(
  model = m1.fr.T1.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=fr, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  ggtitle("Day-10") +
  ylim(0.25,1.75)+
  ylab(expression("Freshness"~"("~italic(alpha)~":"~italic(beta)~")"))+
  xlab("") +
  Fig.formatting

fr.T1.mod.plot

########## T2
m1.fr.T2.gam=gam(fr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")

m2.fr.T2.gam=gam(fr ~  Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")

m3.fr.T2.gam=gam(fr ~  s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")

T2.fr.AIC<-AIC(m1.fr.T2.gam, m2.fr.T2.gam, m3.fr.T2.gam)
T2.fr.AIC
##                    df        AIC
## m1.fr.T2.gam 6.756227 -102.78590
## m2.fr.T2.gam 5.460197 -101.45652
## m3.fr.T2.gam 4.296831  -96.00129
summary(m1.fr.T2.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## fr ~ Treatment + s(plant.mass..g, by = Treatment)
## 
## Parametric coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.790989   0.009899  79.907  < 2e-16 ***
## Treatmentunburned -0.039854   0.013999  -2.847  0.00874 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                      edf Ref.df      F p-value   
## s(plant.mass..g):Treatmentburned   2.237  2.756  3.361 0.03249 * 
## s(plant.mass..g):Treatmentunburned 1.000  1.000 13.326 0.00121 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.493   Deviance explained = 56.7%
## -REML = -41.46  Scale est. = 0.0014698  n = 30
T2.fr.anova=anova(m1.fr.T2.gam)
gam.check(m1.fr.T2.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 11 iterations.
## Gradient range [-1.292444e-05,6.966117e-06]
## (score -41.46041 & scale 0.001469797).
## Hessian positive definite, eigenvalue range [1.292412e-05,13.03067].
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                      k'  edf k-index p-value
## s(plant.mass..g):Treatmentburned   9.00 2.24    1.06    0.48
## s(plant.mass..g):Treatmentunburned 9.00 1.00    1.06    0.57
draw(m1.fr.T2.gam)

concrvity(m1.fr.T2.gam)
## # A tibble: 9 × 3
##   type     term                               concurvity
##   <chr>    <chr>                                   <dbl>
## 1 worst    para                                 5   e- 1
## 2 worst    s(plant.mass..g):Treatmentburned     4.60e-25
## 3 worst    s(plant.mass..g):Treatmentunburned   4.58e-25
## 4 observed para                                 5   e- 1
## 5 observed s(plant.mass..g):Treatmentburned     8.83e-30
## 6 observed s(plant.mass..g):Treatmentunburned   9.84e-31
## 7 estimate para                                 5   e- 1
## 8 estimate s(plant.mass..g):Treatmentburned     1.25e-27
## 9 estimate s(plant.mass..g):Treatmentunburned   1.24e-27
par(mfrow = c(2, 2))
plot(m1.fr.T2.gam, all.terms = TRUE, page=1)

# model predictions
fr.diff.T2<-plot_difference(
  m1.fr.T2.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-1,0.2))+
  Fig.formatting

###########  
#plot for the model output on rawdata
fr.T2.mod.plot<-
  plot_smooths(
  model = m1.fr.T2.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=fr, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
 ylim(0.25,1.75)+
  ggtitle("Day-31") +
  ylab("") +
  xlab("plant material (g)") +
  Fig.formatting

fr.T2.mod.plot

########## T3
m1.fr.T3.gam=gam(fr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")

m2.fr.T3.gam=gam(fr ~  Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")

m3.fr.T3.gam=gam(fr ~  s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")

m4.fr.T3.gam=gam(fr ~  Treatment + s(plant.mass..g, m=c(2,0)), subset = Time.point=="T3", data = DOC.df, method = "REML")

T3.fr.AIC<-AIC(m1.fr.T3.gam, m2.fr.T3.gam, m3.fr.T3.gam)
T3.fr.AIC
##                    df        AIC
## m1.fr.T3.gam 8.337961 -104.51469
## m2.fr.T3.gam 6.186728 -108.46478
## m3.fr.T3.gam 4.617005  -86.80427
summary(m2.fr.T3.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## fr ~ Treatment + s(plant.mass..g)
## 
## Parametric coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.815763   0.009067  89.969  < 2e-16 ***
## Treatmentunburned -0.071997   0.012823  -5.615 7.29e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df     F  p-value    
## s(plant.mass..g) 2.631  3.231 8.329 0.000457 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.661   Deviance explained = 70.3%
## -REML = -46.175  Scale est. = 0.0012332  n = 30
T3.fr.anova=anova(m1.fr.T3.gam)
gam.check(m2.fr.T3.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-2.323279e-09,8.61311e-11]
## (score -46.1752 & scale 0.001233207).
## Hessian positive definite, eigenvalue range [0.6071085,13.5516].
## Model rank =  11 / 11 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                    k'  edf k-index p-value
## s(plant.mass..g) 9.00 2.63    1.26    0.91
draw(m2.fr.T3.gam)

concrvity(m2.fr.T3.gam)
## # A tibble: 6 × 3
##   type     term             concurvity
##   <chr>    <chr>                 <dbl>
## 1 worst    para               5   e- 1
## 2 worst    s(plant.mass..g)   4.61e-25
## 3 observed para               5   e- 1
## 4 observed s(plant.mass..g)   2.12e-31
## 5 estimate para               5   e- 1
## 6 estimate s(plant.mass..g)   1.24e-27
par(mfrow = c(2, 2))
plot(m2.fr.T3.gam, all.terms = TRUE, page=1)

# model predictions
fr.diff.T3<-plot_difference(
  m1.fr.T3.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-1,0.2))+
  Fig.formatting

###########  
#plot for the model output on rawdata
fr.T3.mod.plot<-
  plot_smooths(
  model = m2.fr.T3.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),], 
             aes(x=plant.mass..g, y=fr, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
 ylim(0.25,1.75)+
  geom_line(aes(linetype=Treatment))+
  ggtitle("Day-59") +
  ylab("") +
  xlab("") +
  Fig.formatting

########## T4
m1.fr.T4.gam=gam(fr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")

m2.fr.T4.gam=gam(fr ~  Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")

m3.fr.T4.gam=gam(fr ~  s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")

m4.fr.T4.gam=gam(fr ~  Treatment + s(plant.mass..g, m=c(2,0)), subset = Time.point=="T4", data = DOC.df, method = "REML")

T4.fr.AIC<-AIC(m1.fr.T4.gam, m2.fr.T4.gam, m3.fr.T4.gam)
T4.fr.AIC
##                    df        AIC
## m1.fr.T4.gam 9.278312 -103.54444
## m2.fr.T4.gam 7.825392 -104.36395
## m3.fr.T4.gam 3.470877  -93.80125
summary(m2.fr.T4.gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## fr ~ Treatment + s(plant.mass..g)
## 
## Parametric coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.814924   0.009446   86.27  < 2e-16 ***
## Treatmentunburned -0.043554   0.013359   -3.26  0.00331 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                    edf Ref.df     F  p-value    
## s(plant.mass..g) 3.975  4.825 7.147 0.000483 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.593   Deviance explained = 66.3%
## -REML = -42.963  Scale est. = 0.0013384  n = 30
T4.fr.anova=anova(m1.fr.T4.gam)
gam.check(m2.fr.T4.gam, rep=1000)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-3.584282e-06,1.105561e-06]
## (score -42.96326 & scale 0.001338422).
## Hessian positive definite, eigenvalue range [0.2725441,13.66728].
## Model rank =  11 / 11 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                    k'  edf k-index p-value
## s(plant.mass..g) 9.00 3.98    1.08    0.56
draw(m2.fr.T4.gam)

concrvity(m2.fr.T4.gam)
## # A tibble: 6 × 3
##   type     term             concurvity
##   <chr>    <chr>                 <dbl>
## 1 worst    para               5   e- 1
## 2 worst    s(plant.mass..g)   4.61e-25
## 3 observed para               5   e- 1
## 4 observed s(plant.mass..g)   4.41e-31
## 5 estimate para               5   e- 1
## 6 estimate s(plant.mass..g)   1.24e-27
par(mfrow = c(2, 2))
plot(m2.fr.T4.gam, all.terms = TRUE, page=1)

# model predictions
fr.diff.T4<-plot_difference(
  m1.fr.T4.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-1,0.2))+
  Fig.formatting

###########  
#plot for the model output on rawdata
fr.T4.mod.plot<-
  plot_smooths(
  model = m2.fr.T4.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),], 
             aes(x=plant.mass..g, y=fr, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
 ylim(0.25,1.75)+
  ggtitle("Day-89") +
  geom_line(aes(linetype=Treatment))+
  ylab("") +
  xlab("") +
  Fig.formatting

Generate dataframe for model outputs

####Parametric output
fr.para.sums=as.data.frame(rbind(T1.fr.anova$pTerms.table, T2.fr.anova$pTerms.table, T3.fr.anova$pTerms.table, T4.fr.anova$pTerms.table))
names(fr.para.sums)[1]="df/edf"
#create empty column for ref.df
fr.para.sums <- fr.para.sums %>%
  # Creating an empty column:
  add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
  add_column(Ref.df = "-", .after="df/edf")
  

#####Smooth output
fr.smooth.sums=as.data.frame(rbind(T1.fr.anova$s.table, T2.fr.anova$s.table, T3.fr.anova$s.table, T4.fr.anova$s.table))
names(fr.smooth.sums)[1]="df/edf"
#create empty column for ref.df
fr.smooth.sums <- fr.smooth.sums %>%
  # Creating an empty column:
  add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")

#merge the dataframes
fr.mod.allsums=rbind(fr.para.sums,fr.smooth.sums)


#create new ID column
fr.mod.allsums <- fr.mod.allsums %>%
  # Creating an empty column:
  add_column(Variable= c("FR"), .before="Model")

#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/fr.mod.allsums.pdf", width = 8.3, height = 3.8)       # Export PDF
grid.table(fr.mod.allsums, rows= NULL)
dev.off()

mod.fr.rep<-rep(c(
              "Day 10: Treatment with by-factor smooth", 
              "Day 10: Treatment with global smooth", 
              "Day 10: global smooth",
              "Day 31: Treatment with by-factor smooth", 
              "Day 31: Treatment with global smooth", 
              "Day 31: global smooth",
              "Day 59: Treatment with by-factor smooth", 
              "Day 59: Treatment with global smooth", 
              "Day 59: global smooth",
              "Day 89: Treatment with by-factor smooth", 
              "Day 89: Treatment with global smooth", 
              "Day 89: global smooth"))

mod.fr.column<- data.frame(mod.fr.rep)
AIC.fr<-bind_rows(T1.fr.AIC, T2.fr.AIC, T3.fr.AIC, T4.fr.AIC)
AIC.fr.mod<-as.data.frame(cbind(mod.fr.column, AIC.fr))

names(AIC.fr.mod)[1]="Model"

AIC.fr.mod <- AIC.fr.mod %>%
  # Creating an empty column:
  add_column(Variable= c("FR"), .before="Model")

write.csv(AIC.fr.mod, "output/AIC models/AIC.fr.csv")

#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/fr.mod.AIC.pdf", width = 6.3, height = 3.8)       # Export PDF
grid.table(AIC.fr.mod, rows= NULL)
dev.off()

Freshness: compile raw plots and model-diff plots for final figures

fr.alltimes<-plot_grid(
  fr.T1.mod.plot+ theme(legend.position = "none"),
    fr.T2.mod.plot+ theme(legend.position = "none"),
   fr.T3.mod.plot+ theme(legend.position = "none"),
    fr.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
  rel_widths = c(8,8,8,8,3), ncol=5)

fr.alltimes

ggsave("figures/fr.alltimes.mod.pdf", height=4, width=12)

fr.mod.diffs<-plot_grid(
  fr.diff.T1+ theme(legend.position = "none")+ ggtitle("FR - Day-10"),
  fr.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
  fr.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
  fr.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
  rel_widths = c(8,8,8,8), ncol=4)

fr.mod.diffs

ggsave("figures/fr.mod.diffs.pdf", height=4, width=12)

Fluorescence analysis and plots

library(mgcv)
library(tidymv)
######## T0 model
m1.fl.T0.gam=gam(fl ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")

m2.fl.T0.gam=gam(fl ~  Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")

m3.fl.T0.gam=gam(fl ~  s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")

T0.fl.AIC<-AIC(m1.fl.T0.gam, m2.fl.T0.gam, m3.fl.T0.gam)
T0.fl.AIC

summary(m1.fl.T0.gam)
anova.gam(m1.fl.T0.gam)
gam.check(m1.fl.T0.gam, rep=1000)
draw(m1.fl.T0.gam)
concrvity(m1.fl.T0.gam)
par(mfrow = c(2, 2))
plot(m1.fl.T0.gam, all.terms = TRUE, page=1)

# model predictions
fl.diff.T0<-plot_difference(
  m1.fl.T0.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-1.5,0.5))+
  Fig.formatting

###########  
fl.T0.mod.plot<-
  plot_smooths(
  model = m1.fl.T0.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),], 
             aes(x=plant.mass..g, y=fl, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  ggtitle("Day-0") +
  coord_cartesian(ylim = c(0.5,2.5))+
    ylab("FI")+
  xlab("plant material (g)") +
  Fig.formatting

fl.T0.mod.plot

######## T1 model
m1.fl.T1.gam=gam(fl ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")

m2.fl.T1.gam=gam(fl ~  Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")

m3.fl.T1.gam=gam(fl ~  s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")

T1.fl.AIC<-AIC(m1.fl.T1.gam, m2.fl.T1.gam, m3.fl.T1.gam)
T1.fl.AIC

summary(m1.fl.T1.gam)
T1.fl.anova=anova(m1.fl.T1.gam)
T1.fl.anova
gam.check(m1.fl.T1.gam, rep=1000)
draw(m1.fl.T1.gam)
concrvity(m1.fl.T1.gam)
par(mfrow = c(2, 2))
plot(m1.fl.T1.gam, all.terms = TRUE, page=1)

# model predictions
fl.diff.T1<-plot_difference(
  m1.fl.T1.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-1.5,1))+
  Fig.formatting

###########  
#plot for the model output on rawdata
fl.T1.mod.plot<-
  plot_smooths(
  model = m1.fl.T1.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=fl, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) + 
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  ggtitle("Day-10") +
  ylab("FI")+
  xlab("") +
  ylim(0.5,2.5)+
  Fig.formatting

fl.T1.mod.plot

########## T2
m1.fl.T2.gam=gam(fl ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")

m2.fl.T2.gam=gam(fl ~  Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")

m3.fl.T2.gam=gam(fl ~  s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")

T2.fl.AIC<-AIC(m1.fl.T2.gam, m2.fl.T2.gam, m3.fl.T2.gam)
T2.fl.AIC

summary(m1.fl.T2.gam)
T2.fl.anova=anova(m1.fl.T2.gam)
gam.check(m1.fl.T2.gam, rep=1000)
draw(m1.fl.T2.gam)
concrvity(m1.fl.T2.gam)
par(mfrow = c(2, 2))
plot(m1.fl.T2.gam, all.terms = TRUE, page=1)

# model predictions
fl.diff.T2<-plot_difference(
  m1.fl.T2.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-1.5,1))+
  Fig.formatting

###########  
#plot for the model output on rawdata
fl.T2.mod.plot<-
  plot_smooths(
  model = m1.fl.T2.gam,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),], 
             aes(x=plant.mass..g, y=fl, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  ggtitle("Day-31") +
  ylab("") +
  xlab("plant material (g)") +
  coord_cartesian(ylim = c(0.5,2.5))+
  Fig.formatting

fl.T2.mod.plot

########## T3
m1.fl.T3.gam=gam(fl ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")

m2.fl.T3.gam=gam(fl ~  Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")

m3.fl.T3.gam=gam(fl ~  s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")

T3.fl.AIC<-AIC(m1.fl.T3.gam, m2.fl.T3.gam, m3.fl.T3.gam)
T3.fl.AIC

summary(m3.fl.T3.gam)
T3.fl.anova=anova(m1.fl.T3.gam)
gam.check(m3.fl.T3.gam, rep=1000)
draw(m3.fl.T3.gam)
concrvity(m3.fl.T3.gam)
par(mfrow = c(2, 2))
plot(m3.fl.T3.gam, all.terms = TRUE, page=1)

# model predictions
fl.diff.T3<-plot_difference(
  m1.fl.T3.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-1.5,1))+
  Fig.formatting

###########  
#plot for the model output on rawdata
fl.T3.mod.plot<-
  plot_smooths(
  model = m3.fl.T3.gam,
  series = plant.mass..g,
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),], 
             aes(x=plant.mass..g, y=fl, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0.5,2.5))+
  ggtitle("Day-59") +
  ylab("") +
  xlab("") +
  Fig.formatting

########## T4
m1.fl.T4.gam=gam(fl ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")

m2.fl.T4.gam=gam(fl ~  Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")

m3.fl.T4.gam=gam(fl ~  s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")

T4.fl.AIC<-AIC(m1.fl.T4.gam, m2.fl.T4.gam, m3.fl.T4.gam)
T4.fl.AIC

summary(m3.fl.T4.gam)
T4.fl.anova=anova(m1.fl.T4.gam)
gam.check(m3.fl.T4.gam, rep=1000)
draw(m3.fl.T4.gam)
concrvity(m3.fl.T4.gam)
par(mfrow = c(2, 2))
plot(m3.fl.T4.gam, all.terms = TRUE, page=1)

# model predictions
fl.diff.T4<-plot_difference(
  m1.fl.T4.gam,
  series = plant.mass..g,
  difference = list(Treatment = c("burned", "unburned"))
)+
  coord_cartesian(ylim = c(-1.5,1))+
  Fig.formatting

###########  
#plot for the model output on rawdata
fl.T4.mod.plot<-
  plot_smooths(
  model = m3.fl.T4.gam,
  series = plant.mass..g,
) + 
  geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),], 
             aes(x=plant.mass..g, y=fl, color=Treatment)) +
  scale_color_manual(values = c("brown1", "mediumseagreen")) +
    scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0.5,2.5))+
  ggtitle("Day-89") +
  ylab("") +
  xlab("") +
  Fig.formatting

mod.fl.rep<-rep(c(
              "Day 10: Treatment with by-factor smooth", 
              "Day 10: Treatment with global smooth", 
              "Day 10: global smooth",
              "Day 31: Treatment with by-factor smooth", 
              "Day 31: Treatment with global smooth", 
              "Day 31: global smooth",
              "Day 59: Treatment with by-factor smooth", 
              "Day 59: Treatment with global smooth", 
              "Day 59: global smooth",
              "Day 89: Treatment with by-factor smooth", 
              "Day 89: Treatment with global smooth", 
              "Day 89: global smooth"))

mod.fl.column<- data.frame(mod.fl.rep)
AIC.fl<-bind_rows(T1.fl.AIC, T2.fl.AIC, T3.fl.AIC, T4.fl.AIC)
AIC.fl.mod<-as.data.frame(cbind(mod.fl.column, AIC.fl))
AIC.fl.mod
names(AIC.fl.mod)[1]="Model"

AIC.fl.mod <- AIC.fl.mod %>%
  # Creating an empty column:
  add_column(Variable= c("FI"), .before="Model")


write.csv(AIC.fl.mod, "output/AIC models/AIC.fl.csv")

#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/fl.mod.AIC.pdf", width = 6, height = 3.8)       # Export PDF
grid.table(AIC.fl.mod, rows=NULL)
dev.off()

Generate dataframe for model outputs

####Parametric output
fl.para.sums=as.data.frame(rbind(T1.fl.anova$pTerms.table, T2.fl.anova$pTerms.table, T3.fl.anova$pTerms.table, T4.fl.anova$pTerms.table))
names(fl.para.sums)[1]="df/edf"

#create empty column for ref.df
fl.para.sums <- fl.para.sums %>%
  # Creating an empty column:
  add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
  add_column(Ref.df = "-", .after="df/edf")
  

#####Smooth output
fl.smooth.sums=as.data.frame(rbind(T1.fl.anova$s.table, T2.fl.anova$s.table, T3.fl.anova$s.table, T4.fl.anova$s.table))
names(fl.smooth.sums)[1]="df/edf"
#create empty column for ref.df
fl.smooth.sums <- fl.smooth.sums %>%
  # Creating an empty column:
  add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")

#merge the dataframes
fl.mod.allsums=rbind(fl.para.sums,fl.smooth.sums)
rownames(fl.mod.allsums)<-c(1:12)

#create new ID column
fl.mod.allsums <- fl.mod.allsums %>%
  # Creating an empty column:
  add_column(Variable= c("FI"), .before="Model")


pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/fl.mod.allsums.pdf", width = 8.3, height = 3.8)       # Export PDF
grid.table(fl.mod.allsums, rows=NULL)
dev.off()

Fluorescence: compile raw plots and model-diff plots for final figures


fl.alltimes<-plot_grid(
  fl.T1.mod.plot+ theme(legend.position = "none"),
    fl.T2.mod.plot+ theme(legend.position = "none"),
   fl.T3.mod.plot+ theme(legend.position = "none"),
    fl.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
  rel_widths = c(8,8,8,8,3), ncol=5)

fl.alltimes
ggsave("figures/fr.alltimes.mod.pdf", height=4, width=12)

fl.mod.diffs<-plot_grid(
  fl.diff.T1+ theme(legend.position = "none")+ ggtitle("FI - Day-10"),
  fl.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
  fl.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
  fl.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
  rel_widths = c(8,8,8,8), ncol=4)

fl.mod.diffs
ggsave("figures/fl.mod.diffs.pdf", height=4, width=12)

Compile all indices into one plot


indices.alltimes= plot_grid(suva.alltimes, sr.alltimes, hix.alltimes, fr.alltimes, fl.alltimes, nrow=5, rel_widths = c(8,8,8,8,8), labels = "AUTO")

ggsave("figures/indices.alltimes.pdf", height=18, width=14)

indices.diff.plot= plot_grid(suva.mod.diffs, sr.mod.diffs, hix.mod.diffs, fr.mod.diffs, fl.mod.diffs, nrow=5, rel_widths = c(8,8,8,8,8), labels = "auto")

indices.diff.plot
ggsave("figures/indices.diff.plot.pdf", height=18, width=14)

#### Day 0 baseline indicies

day0_index.baseline.plot = plot_grid(
  suva.T0.mod.plot+ theme(legend.position = "none"),
  sr.T0.mod.plot+ theme(legend.position = "none"),
  hix.T0.mod.plot+ theme(legend.position = "none"),
  fr.T0.mod.plot+ theme(legend.position = "none"),
  fl.T0.mod.plot+ theme(legend.position = "none"),
  extract.legend,
  rel_widths = c(8,8,8,8,8,3), nrow=6)
day0_index.baseline.plot
ggsave("figures/day0_index.baseline.plot.pdf", height=18, width=14)

Calculate photo, microbe, and total degradation:

#+UV/+Bio=A
#-UV/+Bio=B
#+UV/-Bio=C
#-UV/-Bio (NoUV,Filtered)=D

#This is calculating % decrease ((control-treatment)/control)*100

#photodegradation

photo.df=data.frame(UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Time.point,
                    UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Start.DOC,
                    UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Start.TN,
                    UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Treatment,
                    UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Tank,
                    UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$plant.mass..g,
                    UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$DOC..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L,
                    100*(UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$DOC..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L)/(UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L), 
                     UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Start.TN,
                    100*(UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Start.TN)/((UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Start.TN)),
                     suva.df[suva.df$Time.point=="T1",]$SUVA,
                     suva.df[suva.df$Time.point=="T1",]$Spectral.Slope.Ratio,
                    index.df[index.df$Time.point=="T1",]$Humification.Index,
                    index.df[index.df$Time.point=="T1",]$Freshness.Index,
                    index.df[index.df$Time.point=="T1",]$Fluroscence.Index)

photo.df$effect="UV-only"
names(photo.df)[1]="Time.point"
names(photo.df)[2]="Start.DOC"
names(photo.df)[3]="Start.TN"
names(photo.df)[4]="Treatment"
names(photo.df)[5]="Tank"
names(photo.df)[6]="plant.mass..g"
names(photo.df)[7]="abs.change.DOC"
names(photo.df)[8]="percent.change.DOC"
names(photo.df)[9]="abs.change.TN"
names(photo.df)[10]="percent.change.TN"
names(photo.df)[11]="SUVA"
names(photo.df)[12]="Spectral.slope"
names(photo.df)[13]="HIX"
names(photo.df)[14]="Fr"
names(photo.df)[15]="FI"


photo.df$mass.remaining.DOC=100-abs(photo.df[photo.df$effect=="UV-only",]$percent.change.DOC)
photo.df$mass.remaining.TN=100-abs(photo.df[photo.df$effect=="UV-only",]$percent.change.TN)

photo_t1.df=subset(photo.df, Time.point %in% c("T1"))

photo.df=photo.df%>%
    mutate_if(is.character,as.factor)

photo_t1.df=photo_t1.df%>%
    mutate_if(is.character,as.factor)

#Microbe decomp
microbe.df=data.frame(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Time.point,
                      UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Start.DOC,
                      UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Start.TN,
                      UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Treatment,
                      UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Tank,
                      UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$plant.mass..g,
                    UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$DOC..mg.L- UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L,
                    100*(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$DOC..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L)/(( UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L)),
                     UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Start.TN,
                    100*(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Start.TN)/((UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Start.TN)),
                     suva.df[suva.df$Time.point=="T1",]$SUVA,
                    suva.df[suva.df$Time.point=="T1",]$Spectral.Slope.Ratio,
                    index.df[index.df$Time.point=="T1",]$Humification.Index,
                    index.df[index.df$Time.point=="T1",]$Freshness.Index,
                    index.df[index.df$Time.point=="T1",]$Fluroscence.Index)

microbe.df$effect="microbes-only"
names(microbe.df)[1]="Time.point"
names(microbe.df)[2]="Start.DOC"
names(microbe.df)[3]="Start.TN"
names(microbe.df)[4]="Treatment"
names(microbe.df)[5]="Tank"
names(microbe.df)[6]="plant.mass..g"
names(microbe.df)[7]="abs.change.DOC"
names(microbe.df)[8]="percent.change.DOC"
names(microbe.df)[9]="abs.change.TN"
names(microbe.df)[10]="percent.change.TN"
names(microbe.df)[11]="SUVA"
names(microbe.df)[12]="Spectral.slope"
names(microbe.df)[13]="HIX"
names(microbe.df)[14]="Fr"
names(microbe.df)[15]="FI"
microbe.df$mass.remaining.DOC=100-abs(microbe.df[microbe.df$effect=="microbes-only",]$percent.change.DOC)
microbe.df$mass.remaining.TN=100-abs(microbe.df[microbe.df$effect=="microbes-only",]$percent.change.TN)

microbe_t1.df=subset(microbe.df, Time.point %in% c("T1"))

microbe.df=microbe.df%>%
    mutate_if(is.character,as.factor)

microbe_t1.df=microbe_t1.df%>%
    mutate_if(is.character,as.factor)

#Total decomp
total.df=data.frame(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Time.point,
                    UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Start.DOC,
                    UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Start.TN,
                    UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Treatment,
                    UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Tank,
                    UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$plant.mass..g,
                    UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$DOC..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L,
                     100*(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$DOC..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L)/(abs(UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L)),
                      UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Start.TN,
                     100*(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Start.TN)/(abs(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Start.TN)),
                     suva.df[suva.df$Time.point=="T1",]$SUVA,
                    suva.df[suva.df$Time.point=="T1",]$Spectral.Slope.Ratio,
                     index.df[index.df$Time.point=="T1",]$Humification.Index,
                    index.df[index.df$Time.point=="T1",]$Freshness.Index,
                    index.df[index.df$Time.point=="T1",]$Fluroscence.Index)
                     
total.df$effect="UV and microbes"
names(total.df)[1]="Time.point"
names(total.df)[2]="Start.DOC"
names(total.df)[3]="Start.TN"
names(total.df)[4]="Treatment"
names(total.df)[5]="Tank"
names(total.df)[6]="plant.mass..g"
names(total.df)[7]="abs.change.DOC"
names(total.df)[8]="percent.change.DOC"
names(total.df)[9]="abs.change.TN"
names(total.df)[10]="percent.change.TN"
names(total.df)[11]="SUVA"
names(total.df)[12]="Spectral.slope"
names(total.df)[13]="HIX"
names(total.df)[14]="Fr"
names(total.df)[15]="FI"
total.df$mass.remaining.DOC=100-abs(total.df[total.df$effect=="UV and microbes",]$percent.change.DOC)
total.df$mass.remaining.TN=100-abs(total.df[total.df$effect=="UV and microbes",]$percent.change.TN)

total_t1.df=subset(total.df, Time.point %in% c("T1"))

total.df=total.df%>%
    mutate_if(is.character,as.factor)

total_t1.df=total_t1.df%>%
    mutate_if(is.character,as.factor)

#merge degra.df time points into one dataframe


degra.df=rbind(photo.df, microbe.df, total.df)
degra.df=pivot_longer(degra.df, cols=c("effect"), names_to="effect", values_to="Type")

degra_t1.df=rbind(photo_t1.df, microbe_t1.df, total_t1.df)
degra_t1.df=pivot_longer(degra_t1.df, cols=c("effect"), names_to="effect", values_to="Type")

#make sure we have factors
degra.df=degra.df%>%
    mutate_if(is.character,as.factor)

degra_t1.df=degra_t1.df%>%
    mutate_if(is.character,as.factor)

#seperate degra.df for burned and unburned tanks
degra_t1.burned.df=degra_t1.df[degra_t1.df$Treatment=="burned",]
degra_t1.unburned.df=degra_t1.df[degra_t1.df$Treatment=="unburned",]

UVBIO models and plots

library(mgcv)

###Just UV
m1.photo.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g, by= Treatment), data=photo.df, method="REML", select= T, subset= Time.point=="T1")

m2.photo.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g), data=photo.df, method="REML", select= T, subset= Time.point=="T1")

m3.photo.perc=gam(percent.change.DOC ~ s(plant.mass..g), data=photo.df, method="REML", select= T, subset= Time.point=="T1")

m1.photo.abs=gam(abs.change.DOC ~ Treatment +s(plant.mass..g, by= Treatment), data=photo.df, method="REML", select= T, subset= Time.point=="T1")

m2.photo.abs=gam(abs.change.DOC ~ Treatment +s(plant.mass..g), data=photo.df, method="REML", select= T, subset= Time.point=="T1")

m3.photo.abs=gam(abs.change.DOC ~ s(plant.mass..g), data=photo.df, method="REML", select= T, subset= Time.point=="T1")
photo.AIC=AIC(m1.photo.perc, m2.photo.perc, m3.photo.perc)
photo.AIC=AIC(m1.photo.abs, m2.photo.abs, m3.photo.abs)


photo.AIC

summary(m1.photo.perc)
photo.anova=anova(m1.photo.perc)
photo.anova
gam.check(m3.photo.perc, rep=1000)
draw(m3.photo.perc)
concrvity(m3.photo.perc)
par(mfrow = c(2, 2))
plot(m3.photo.perc, all.terms = TRUE, page=1)


#plot for the model output on rawdata
photo.mod.plot<-
  plot_smooths(
  model = m1.photo.perc,
  series = plant.mass..g,
  comparison = Treatment
) + 
  geom_point(data=photo.df[(photo.df$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=percent.change.DOC, color=Treatment)) +
   scale_color_manual(values = c("brown1", "mediumseagreen")) +
  geom_hline(yintercept=0, linetype="longdash", color = "gray")+
  ggtitle("UV-only") +
  ylab("Change DOC (%)") +
  xlab("plant material (g)") +
    coord_cartesian(ylim = c(-30,10))+
  Fig.formatting

photo.mod.plot

###Just microbes
m1.microbe.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g, by= Treatment), data=microbe.df, method="REML", select= T, subset= Time.point=="T1")

m2.microbe.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g), data=microbe.df, method="REML", select= T, subset= Time.point=="T1")

m3.microbe.perc=gam(percent.change.DOC ~ s(plant.mass..g), data=microbe.df, method="REML", select= T, subset= Time.point=="T1")

microbe.AIC=AIC(m1.microbe.perc, m2.microbe.perc, m3.microbe.perc)
microbe.AIC

summary(m3.microbe.perc)
microbe.anova=anova(m1.microbe.perc)
gam.check(m3.microbe.perc, rep=1000)
draw(m3.microbe.perc)
concrvity(m3.microbe.perc)
par(mfrow = c(2, 2))
plot(m3.microbe.perc, all.terms = TRUE, page=1)

#plot for the model output on rawdata
microbe.mod.plot<-
  plot_smooths(
  model = m3.microbe.perc,
  series = plant.mass..g
) + 
  geom_point(data=microbe.df[(microbe.df$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=percent.change.DOC, color=Treatment)) +
   scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ggtitle("microbes-only") +
   geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab("Change DOC (%)") +
   coord_cartesian(ylim = c(-30,10))+
  xlab("plant material (g)") +
  Fig.formatting

 coord_cartesian(ylim = c(-30,0))


###Just total decomp
m1.total.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g, by= Treatment), data=total.df, method="REML", select= T, subset= Time.point=="T1")

m2.total.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g), data=total.df, method="REML", select= T, subset= Time.point=="T1")

m3.total.perc=gam(percent.change.DOC ~ s(plant.mass..g), data=total.df, method="REML", select= T, subset= Time.point=="T1")

m4.total.perc=gam(percent.change.DOC ~ s(plant.mass..g, m=c(2,0)), data=total.df, method="REML", select= T, subset= Time.point=="T1")

total.AIC=AIC(m1.total.perc, m2.total.perc, m3.total.perc)
total.AIC

summary(m3.total.perc)
total.anova=anova(m1.total.perc)
gam.check(m3.total.perc, rep=1000)
draw(m3.total.perc)
concrvity(m3.total.perc)
par(mfrow = c(2, 2))
plot(m3.total.perc, all.terms = TRUE, page=1)

#plot for the model output on rawdata
total.mod.plot<-
  plot_smooths(
  model = m3.total.perc,
  series = plant.mass..g
) + 
  geom_point(data=total.df[(total.df$Time.point=="T1"),], 
             aes(x=plant.mass..g, y=percent.change.DOC, color=Treatment)) +
   scale_color_manual(values = c("brown1", "mediumseagreen")) +
  ggtitle("UV and microbes") +
   geom_hline(yintercept=0, linetype="longdash", color = "gray") +
  ylab("Change DOC (%)") +
  coord_cartesian(ylim = c(-30,10))+
  xlab("plant material (g)") +
  Fig.formatting

##combine AIC
mod.degra.rep<-rep(c("UV-only: Treatment with by-factor smooth", 
              "UV-only: Treatment with global smooth", 
              "UV-only: global smooth",
              "Microbes-only: Treatment with by-factor smooth", 
              "Microbes-only: Treatment with global smooth", 
              "Microbes-only: global smooth",
              "UV and microbes: Treatment with by-factor smooth", 
              "UV and microbes: Treatment with global smooth", 
              "UV and microbes: global smooth"))

mod.degra.column<- data.frame(mod.degra.rep)
AIC.degra<-bind_rows(photo.AIC, microbe.AIC, total.AIC)
AIC.degra.mod<-as.data.frame(cbind(mod.degra.column, AIC.degra))
AIC.degra.mod
names(AIC.degra.mod)[1]="Model"

AIC.degra.mod <- AIC.degra.mod %>%
  # Creating an empty column:
  add_column(Variable= c("-"), .before="Model")


write.csv(AIC.degra.mod, "output/AIC models/AIC.degra.csv")

#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/degra.mod.AIC.pdf", width = 6.5, height = 2.8)       # Export PDF
grid.table(AIC.degra.mod, rows=NULL)
dev.off()

calculate effect size for degradation pathways

library(effectsize)
r_squared_microbes <- summary(m1.microbe.perc)$r.sq
r_squared_total <- summary(m1.total.perc)$r.sq
r_squared_microbes
## [1] 0.3219554
r_squared_total
## [1] 0.4629789

Generate dataframe for model outputs

####Parametric output
degra.para.sums=as.data.frame(rbind(photo.anova$pTerms.table, microbe.anova$pTerms.table, total.anova$pTerms.table))
names(degra.para.sums)[1]="df/edf"

#create empty column for ref.df
degra.para.sums <- degra.para.sums %>%
  # Creating an empty column:
  add_column(Model= c("UV-only Model: Treatment","Microbes-only Model: Treatment","UV and microbes Model: Treatment"), .before="df/edf")%>%
  add_column(Ref.df = "-", .after="df/edf")
  

#####Smooth output
degra.smooth.sums=as.data.frame(rbind(photo.anova$s.table, microbe.anova$s.table, total.anova$s.table))
names(degra.smooth.sums)[1]="df/edf"

#create empty column for ref.df
degra.smooth.sums <- degra.smooth.sums %>%
  # Creating an empty column:
  add_column(Model= c("UV-only Model: burned smooth","UV-only Model: unburned smooth","Microbes-only Model: burned smooth","Microbes-only Model: unburned smooth","UV and microbes Model: burned smooth","UV and microbes Model: unburned smooth"), .before="df/edf")

#merge the dataframes
degra.mod.allsums=rbind(degra.para.sums,degra.smooth.sums)

degra.mod.allsums <- degra.mod.allsums %>%
  # Creating an empty column:
  add_column(Variable= c("-"), .before="Model")


write.csv(degra.mod.allsums,"/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/degra.mod.allsums.csv")

#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/degra.mod.allsums.pdf", width = 8.5, height = 2.8)       # Export PDF
grid.table(degra.mod.allsums, rows=NULL)
dev.off()

degra mod plots

degra.mod.plots = plot_grid(
  photo.mod.plot + theme(legend.position = "none"),
  microbe.mod.plot + theme(legend.position = "none"),
  total.mod.plot + theme(legend.position = "none" ),
  extract.legend,
  rel_widths = c(8,8,8,3), nrow = 1)

degra.mod.plots

ggsave("figures/degra.mod.plots.pdf", height=4, width=10)

Sage and Willow Decomposition Analysis

#load file
plant.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania2/data/decomp.weights.csv")
plant.df=plant.df[-c(61:82), ]
plant.df=plant.df[-c(6, 8,9)]

#make sure we have factors
plant.df=plant.df%>%
    mutate_if(is.character,as.factor)

#add propotion column so it can be arcsine transformed
plant.df$perc.loss=(100*(plant.df[plant.df$plant.type=="sage" | plant.df$plant.type=="willow",]$plant.specific.mass..g-plant.df[plant.df$plant.type=="sage" | plant.df$plant.type=="willow",]$dry.mass..g.bag.corrected)/(plant.df[plant.df$plant.type=="sage" | plant.df$plant.type=="willow",]$plant.specific.mass..g))

#remove NaN
plant.df[is.na(plant.df)] <- 0 

#calculate mass remaining
plant.df$mass.remaining=100-plant.df$perc.loss
hist(plant.df$mass.remaining)

#add delta
plant.df$delta.loss= plant.df$perc.loss/100
  
#arcsine
plant.df$delta.trans=asin(sqrt(plant.df$delta.loss))
hist(plant.df$delta.trans)


#exclude control tank outlier

plant.df2=data.frame(plant.df[plant.df$plant.specific.mass..g > 0,])


###Sage
m1.sage.gam=gam(delta.trans ~ treatment +s(plant.specific.mass..g, by= treatment), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")



m2.sage.gam=gam(delta.trans ~ treatment +s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")

m3.sage.gam=gam(delta.trans ~ s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")

m4.sage.gam=gam(perc.loss ~ treatment +s(plant.specific.mass..g, by= treatment), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")

hist(m3.sage.gam$residuals)

m5.sage.gam=gam(perc.loss ~ treatment +s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")

m6.sage.gam=gam(perc.loss ~ s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")

sage.AIC=AIC(m4.sage.gam, m5.sage.gam, m6.sage.gam)
sage.AIC

sage.anova=anova(m4.sage.gam)
sage.anova
gam.check(m6.sage.gam, rep=1000)

mean(plant.df[plant.df$plant.type=="sage",]$delta.trans)
# global smooth

# model predictions
sage.diff.mod<-plot_difference(
  m1.sage.gam,
  series = plant.specific.mass..g,
  difference = list(treatment = c("burned", "unburned"))
)

sage.diff.mod

summary(m4.sage.gam)
anova.gam(m3.sage.gam)
gam.check(m4.sage.gam, rep=1000)
draw(m3.sage.gam)
concrvity(m3.sage.gam)
par(mfrow = c(2, 2))
plot(m3.sage.gam, all.terms = TRUE, page=1)
qq_plot(m3.sage.gam)

#plot for the model output on rawdata


sage.mod.plot<-
  plot_smooths(
  model = m4.sage.gam,
  series = plant.specific.mass..g,
  comparison= treatment
) + 
  geom_point(data=plant.df2[(plant.df2$plant.type=="sage"),], 
             aes(x=plant.specific.mass..g, y=perc.loss, color=treatment)) +
   scale_color_manual(values = c("brown1", "mediumseagreen")) +
  scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  ggtitle("Sage") +
 ylab("Mass loss (%)") +
    coord_cartesian(ylim = c(0,80))+
  xlab("Plant material (g)") +
  Fig.formatting

sage.mod.plot


### Willow
m1.willow.gam=gam(delta.trans ~ treatment +s(plant.specific.mass..g, by= treatment), data=plant.df2, method="REML", select= T, subset= plant.type=="willow")


m2.willow.gam=gam(delta.trans ~ treatment +s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="willow")

m3.willow.gam=gam(delta.trans ~ s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="willow")

m4.willow.gam=gam(perc.loss ~ treatment +s(plant.specific.mass..g, by= treatment), data=plant.df2, method="REML", select= T, subset= plant.type=="willow")


m5.willow.gam=gam(perc.loss ~ treatment +s(plant.specific.mass..g), data=plant.df, method="REML", select= T, subset= plant.type=="willow")

m6.willow.gam=gam(perc.loss ~ s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="willow")

willow.AIC=AIC( m4.willow.gam,m5.willow.gam,m6.willow.gam)
willow.AIC
# treatment with global smooth

summary(m4.willow.gam)
willow.anova=anova(m4.willow.gam)
willow.anova
gam.check(m3.willow.gam, rep=1000)
draw(m2.willow.gam)
concrvity(m2.willow.gam)
par(mfrow = c(2, 2))
plot(m2.willow.gam, all.terms = TRUE, page=1)

# model predictions
willow.diff.mod<-plot_difference(
  m1.willow.gam,
  series = plant.specific.mass..g,
  difference = list(treatment = c("burned", "unburned"))
)
willow.diff.mod

#plot for the model output on rawdata
willow.mod.plot<-
  plot_smooths(
  model = m4.willow.gam,
  series = plant.specific.mass..g,
  comparison = treatment
) + 
  geom_point(data=plant.df2[(plant.df2$plant.type=="willow"),], 
             aes(x=plant.specific.mass..g, y=perc.loss, color=treatment)) +
   scale_color_manual(values = c("brown1", "mediumseagreen")) +
     scale_fill_manual(values = c("brown1", "mediumseagreen")) +
  ggtitle("Willow") +
  ylab("Mass loss (%)") +
  xlab("Plant material (g)") +
  coord_cartesian(ylim = c(0,80))+
  Fig.formatting


willow.mod.plot


##combine AIC

mod.decomp.rep<-rep(c("Sage: Treatment with by-factor smooth", 
              "Sage: Treatment with global smooth", 
              "Sage: global smooth",
              "Willow: Treatment with by-factor smooth", 
              "Willow: Treatment with global smooth", 
              "Willow: global smooth"))

mod.decomp.column<- data.frame(mod.decomp.rep)
AIC.decomp<-bind_rows(sage.AIC, willow.AIC)
AIC.decomp.mod<-as.data.frame(cbind(mod.decomp.column, AIC.decomp))

names(AIC.decomp.mod)[1]="Model"


write.csv(AIC.decomp.mod, "output/AIC models/AIC.decomp.csv")

#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/decomp.mod.AIC.pdf", width = 5, height = 2)       # Export PDF
grid.table(AIC.decomp.mod, rows=NULL)
dev.off()

Generate dataframe for model outputs

####Parametric output
decomp.para.sums=as.data.frame(rbind(sage.anova$pTerms.table, willow.anova$pTerms.table))
names(decomp.para.sums)[1]="df/edf"

#create empty column for ref.df
decomp.para.sums <- decomp.para.sums %>%
  # Creating an empty column:
  add_column(Model= c("Sage Model: Treatment","Willow Model: Treatment"), .before="df/edf")%>%
  add_column(Ref.df = "-", .after="df/edf")
  

#####Smooth output
decomp.smooth.sums=as.data.frame(rbind(sage.anova$s.table, willow.anova$s.table))
names(decomp.smooth.sums)[1]="df/edf"

#create empty column for ref.df
decomp.smooth.sums <- decomp.smooth.sums %>%
  # Creating an empty column:
  add_column(Model= c("Sage Model: burned smooth","Sage Model: unburned smooth","Willow Model: burned smooth","Willow Model: unburned smooth"), .before="df/edf")

#merge the dataframes
decomp.mod.allsums=rbind(decomp.para.sums,decomp.smooth.sums)



#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/decomp.mod.allsums.pdf", width = 7, height = 2)       # Export PDF
grid.table(decomp.mod.allsums, rows=NULL)
dev.off()

comparing sage and willow by fire treatment

#compare mass loss between sage and willow by fire treatment---burned
m1.compare.burned.gam=gam(perc.loss ~ plant.type + s(plant.specific.mass..g, by=plant.type), data=plant.df2, methods= "REML", select= T, subset= treatment == "burned")

m2.compare.burned.gam=gam(perc.loss ~ plant.type + s(plant.specific.mass..g), data=plant.df2, methods= "REML", select= T, subset= treatment == "burned")

m3.compare.burned.gam=gam(perc.loss ~ s(plant.specific.mass..g), data=plant.df2, methods= "REML", select= T, subset= treatment == "burned")

compare.burned.AIC=AIC(m1.compare.burned.gam, m2.compare.burned.gam, m3.compare.burned.gam)
compare.burned.anova=anova(m1.compare.burned.gam)
compare.burned.anova
gam.check(m1.compare.burned.gam, rep=1000)

# model predictions
compare.burned.diff.mod<-plot_difference(
  m1.compare.burned.gam,
  series = plant.specific.mass..g,
  difference = list(plant.type = c("sage", "willow"))
)
compare.burned.diff.mod

#plot for the model output on rawdata
compare.burned.mod.plot<-
  plot_smooths(
  model = m2.compare.burned.gam,
  series = plant.specific.mass..g,
  comparison = plant.type
) + 
  geom_point(data=plant.df2[(plant.df2$treatment=="burned"),], 
             aes(x=plant.specific.mass..g, y=perc.loss, color=plant.type)) +
   scale_color_manual(values = c("mediumpurple4", "darkgoldenrod4")) +
     scale_fill_manual(values = c("mediumpurple4", "darkgoldenrod4")) +
  ggtitle("Burned") +
  geom_line(aes(linetype= plant.type))+
  ylab("Mass loss (%)") +
  xlab("Plant material (g)") +
    coord_cartesian(ylim = c(0,80))+
  Fig.formatting

compare.burned.mod.plot


#compare mass loss between sage and willow by fire treatment---unburned
m1.compare.unburned.gam=gam(perc.loss ~ plant.type + s(plant.specific.mass..g, by=plant.type), data=plant.df2, methods= "REML", select= T, subset= treatment == "unburned")

m2.compare.unburned.gam=gam(perc.loss ~ plant.type + s(plant.specific.mass..g), data=plant.df2, methods= "REML", select= T, subset= treatment == "unburned")

m3.compare.unburned.gam=gam(perc.loss ~ s(plant.specific.mass..g), data=plant.df2, methods= "REML", select= T, subset= treatment == "unburned")

compare.unburned.AIC=AIC(m1.compare.unburned.gam, m2.compare.unburned.gam, m3.compare.unburned.gam)
compare.unburned.anova=anova(m1.compare.unburned.gam)
compare.unburned.anova
summary(m3.compare.unburned.gam)
gam.check(m3.compare.unburned.gam, rep=1000)

mean(plant.df2[plant.df2$plant.type=="sage",]$perc.loss)
mean(plant.df2[plant.df2$plant.type=="willow",]$perc.loss)
# model predictions
compare.burned.diff.mod<-plot_difference(
  m1.compare.burned.gam,
  series = plant.specific.mass..g,
  difference = list(plant.type = c("sage", "willow"))
)
compare.burned.diff.mod

#plot for the model output on rawdata
compare.unburned.mod.plot<-
  plot_smooths(
  model = m1.compare.unburned.gam,
  series = plant.specific.mass..g,
  comparison = plant.type
) + 
  geom_point(data=plant.df2[(plant.df2$treatment=="unburned"),], 
             aes(x=plant.specific.mass..g, y=perc.loss, color=plant.type)) +
   scale_color_manual(values = c("mediumpurple4", "darkgoldenrod4")) +
     scale_fill_manual(values = c("mediumpurple4", "darkgoldenrod4")) +
  ggtitle("Unburned") +
  ylab("Mass Loss (%)") +
  xlab("Plant material (g)") +
  coord_cartesian(ylim = c(0,80))+
  Fig.formatting

compare.unburned.mod.plot

##combine AIC

mod.comapre.decomp.rep<-rep(c("Burned: plant type with by-factor smooth", 
              "Burned: plant type with global smooth", 
              "Burned: global smooth",
              "Unburned: plant type with by-factor smooth", 
              "Unburned: plant type with global smooth", 
              "Unburned: global smooth"))

mod.compare.decomp.column<- data.frame(mod.comapre.decomp.rep)
AIC.compare.decomp<-bind_rows(compare.burned.AIC, compare.unburned.AIC)
AIC.compare.decomp.mod<-as.data.frame(cbind(mod.compare.decomp.column, AIC.compare.decomp))

names(AIC.compare.decomp.mod)[1]="Model"


write.csv(AIC.compare.decomp.mod, "output/AIC models/AIC.comapre.decomp.csv")

#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/decomp.compare.mod.AIC.pdf", width = 5.3, height = 2)       # Export PDF
grid.table(AIC.compare.decomp.mod, rows=NULL)
dev.off()

compile raw plots and model-diff plots for final figures

decomp.compare.mod.plot<-plot_grid(
compare.burned.mod.plot+ theme(legend.position = "none"),
compare.unburned.mod.plot+ theme(legend.position = "none"),
ncol=2)

decomp.compare.mod.plot

ggsave("figures/decomp.comapre.mod.plot.pdf", height=4, width=12)

Generate dataframe for model outputs

####Parametric output
decomp.compare.para.sums=as.data.frame(rbind(compare.burned.anova$pTerms.table, compare.unburned.anova$pTerms.table))
names(decomp.compare.para.sums)[1]="df/edf"

#create empty column for ref.df
decomp.compare.para.sums <- decomp.compare.para.sums %>%
  # Creating an empty column:
  add_column(Model= c("Burned: Plant type","Unburned: Plant type"), .before="df/edf")%>%
  add_column(Ref.df = "-", .after="df/edf")
  

#####Smooth output
decomp.compare.smooth.sums=as.data.frame(rbind(compare.burned.anova$s.table, compare.unburned.anova$s.table))
names(decomp.compare.smooth.sums)[1]="df/edf"

#create empty column for ref.df
decomp.compare.smooth.sums <- decomp.compare.smooth.sums %>%
  # Creating an empty column:
  add_column(Model= c("Burned: sage smooth","Burned: willow smooth","Unburned Model: sage smooth","Unburned Model: willow smooth"), .before="df/edf")

#merge the dataframes
decomp.compare.mod.allsums=rbind(decomp.compare.para.sums,decomp.compare.smooth.sums)



#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/decomp.compare.mod.allsums.pdf", width = 6.2, height = 2)       # Export PDF
grid.table(decomp.compare.mod.allsums, rows=NULL)
dev.off()

Comnine AIC and Mod sums for decomposition

library(gridExtra)
#merge all the dataframes
all.decomp.mod.sums=rbind(decomp.mod.allsums,decomp.compare.mod.allsums)

#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/all.decomp.mod.sums.pdf", width = 7, height = 3.7)       # Export PDF
grid.table(all.decomp.mod.sums, rows=NULL)
dev.off()

#merge all the dataframes
all.decomp.mod.AIC=rbind(AIC.decomp.mod, AIC.compare.decomp.mod)

#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/all.decomp.mod.AIC.pdf", width = 5.3, height = 3.7)       # Export PDF
grid.table(all.decomp.mod.AIC, rows=NULL)
dev.off()

compile raw plots and model-diff plots for final figures

decomp.mod.plot<-plot_grid(
sage.mod.plot+ theme(legend.position = "none"),
willow.mod.plot+ theme(legend.position = "none"),
  rel_widths = c(8,8), ncol=2)

decomp.mod.plot

ggsave("figures/decomp.mod.plot.pdf", height=4, width=12)

compile raw plots and model-diff plots for final figures


all.decomp.mod.plot<-plot_grid(
decomp.mod.plot+ theme(legend.position = "none"),
decomp.compare.mod.plot+ theme(legend.position = "none"),
nrow=2)

all.decomp.mod.plot

ggsave("figures/all.decomp.mod.plot.pdf", height=8, width=8)